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Abstract 

An excess in 7 -rays emanating from the galactic centre has recently been observed in the Fermi- 
LAT data. This signal can be interpreted as resulting from WIMP annihilation, with the 
spectrum well-fit by dark matter annihilating dominantly into either bottom-quark or Higgs 
pairs. Supersymmetric models provide a well-motivated framework to study the implications of 
this signal in these channels. With a neutralino dark matter candidate, the 7 -ray excess can¬ 
not be easily accommodated in the minimal supersymmetric model, which in any case requires 
tuning below the percent level to explain the observed Higgs mass. Instead we are naturally led 
to consider the next-to-minimal model with a singlet superfield. This not only allows for the 
annihilation channel into bottom-quark pairs to be implemented, but also provides new possi¬ 
bilities for annihilation into Higgs-pseudoscalar pairs. We show that the fit to the 7 -ray excess 
for the Higgs-pseudoscalar channel can be just as good as for annihilation into bottom-quark 
pairs. Moreover, in the parameter range of interest, the next-to-minimal supersymmetric model 
solves the /x-problem and can explain the 125 GeV Higgs mass with improved naturalness. We 
also consider an extension by adding a right-handed neutrino superfield with the right-handed 
sneutrino acting as a dark matter candidate. Interestingly, this allows for the annihilation into 
pseudoscalar pairs which also provide a good fit to the 7 -ray excess. Furthermore, in the case of 
a neutralino LSP, the late decay of a sneutrino NLSP can non-thermally produce the observed 
relic abundance. Finally, the WIMP annihilation into scalar pairs allows for the possibility of 
detecting the Higgs or pseudoscalar decay into two photons, providing a smoking-gun signal of 
the model. 
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1 Introduction 


Cosmic-ray experiments are a promising way to search for dark matter (DM). In particular, the 
satellite-based experiment Fermi-LAT is able to measure the 7 -ray sky with unprecedented preci¬ 
sion. If DM annihilates into photons with energies from 20 MeV to 300 GeV, an imprint can be left 
on these measurements. However, the astrophysical 7 -ray background has to be carefully modelled 
in order to disentangle a putative DM signal. To be consistent with DM, any excess over this 
background needs to have a morphology which is consistent with the range of DM density profiles 
indicated by observations and simulations of structure formation. 

Intriguingly, an excess from the galactic centre with such a morphology has been identified 
in the Fermi-LAT data [1—10]. The spectrum of this Galactic Center Excess (GCE) peaks at 
energies of about 3 GeV. It can be well fit by DM with mass in the range 30 — 70 GeV which 
annihilates dominantly into bb -pairs (or lighter quarks). Interestingly, the annihilation cross section 
that is then required to explain the GCE is of the right size to account for the DM density from 
thermal freeze-out. The systematic uncertainties of the spectrum from modelling the astrophysical 
7 -ray background have recently been studied in Refs. [10; 11]. It was found that, taking the 
estimated uncertainty in the high-energy tail of the spectrum into account, DM annihilating to 
Higgs pairs close to threshold also gives a good fit [11; 12]. Moreover, the Fermi-LAT collaboration 
has recently confirmed the GCE and presented preliminary spectra (although without information 
on the errors) [13]. It was subsequently pointed out [12] that this may also allow for annihilation into 
top-, W- and Z-pairs to reproduce the spectrum well. However, care must be taken in interpreting 
these results as they are still preliminary and little detail of the modelling and uncertainties have 
been publicly released. We will therefore focus on the best-fit channels, to bb- and Higgs-pairs 
(and possibly other scalars), for the spectrum from Ref. [11]. It is then interesting to study well- 
motivated DM models in which these channels arise naturally (see e.g. Refs. [14-21] for earlier 
studies of DM models to explain the GCE). 

A weakly interacting massive particle (WIMP) has the right properties to account for the 
observed DM from thermal freeze-out and is naturally obtained in supersymmetric models with 
P-parity. An excellent DM candidate is the lightest neutralino. The annihilation cross section of 
neutralinos into Higgs-pairs is p -wave suppressed [22] . Due to the small velocities in our galaxy, this 
is therefore not large enough to generate the GCE. We will instead consider the related annihilation 
channel of a Higgs and a Higgs-sector pseudoscalar, which proceeds at the s-wave level (for two 
pseudoscalars in the final state, the cross section is again p -wave suppressed). Since this channel 
has not been considered before (with the exception of Ref. [23] which, however, focused on a singlet¬ 
like scalar instead of the Higgs in the final state), we have performed a fit using the spectra and 
error estimates from Ref. [10]. For comparison, we have also performed fits for the annihilation 
into bb- and Higgs-pairs. We find that an acceptable fit is obtained for pseudoscalar masses up to 
around 150 GeV. In particular, in the region where the pseudoscalar is lighter than the Higgs, 
the fit is improved compared to the annihilation channel into Higgs-pairs. For a sufficiently light 
pseudoscalar, the fit becomes even better than that for the 66 -channel. 

Similarly a light pseudoscalar is required to mediate the 66 -channel. Even though the annihi¬ 
lation of neutralinos into quarks can be mediated by i/u-channel exchange of squarks, the cross 
section for this process is typically not large enough to generate the GCE due to collider constraints 
on the squark masses [24]. In the s-channel, possible mediators are the Z-boson and CP-even and 
CP-odd scalars from the Higgs sector. The Z-boson and CP-even scalars lead to annihilation 
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cross sections that are respectively helicity and p -wave suppressed and can thus not significantly 
contribute to the GCE either. Higgs-sector pseudoscalars will therefore be the dominant mediators 
in the s-channel. As was previously observed [17; 25], the corresponding annihilation cross section 
needs to be resonantly enhanced in order to generate a 7 -ray flux consistent with the GCE. Given 
that the DM mass that best fits the GCE for the fob-channel is in the range 30 — 70GeV [1-10], 
this means that the pseudoscalar mediator needs to be lighter than about 140 GeV. 

The light pseudoscalar required for the Higgs-pseudoscalar and fofo-channels is, however, difficult 
to obtain in the minimal supersymmetric standard model (MSSM ). 1 Indeed, CMS excludes an 
MSSM Higgs sector with a pseudoscalar mass below 180 GeV for values of tan (3 required to obtain 
the observed Higgs mass for stop masses around 1 TeV [27] . Even though this bound may be relaxed 
for heavier stops (see, however, Ref. [28]), this leads to excessive fine-tuning (which is already high 
in the MSSM). In addition, a light pseudoscalar implies that also the second CP-even Higgs and 
the charged Higgs are light. This results in large deviations in the Higgs couplings from standard 
model (SM) values (see e.g. Refs. [29; 30]) and large contributions to flavour-changing processes 
(see e.g. Ref. [31]), respectively, both in conflict with experiment. There may be ways out to avoid 
these constraints but they would require accidental cancellations and thus be likely contrived. 

Given these problems in the MSSM, we will extend the Higgs sector by a singlet superfield 
and consider the next-to-minimal supersymmetric standard model (NMSSM). The lightest pseu¬ 
doscalar is then generically an admixture of the MSSM pseudoscalar and the CP-odd component 
of the singlet. The mass of the lightest pseudoscalar can then be decoupled from the masses of the 
(doublet-like) CP-even and charged Higgs states, allowing the collider constraints in the MSSM 
to be sufficiently relaxed. In addition, the coupling of the Higgs to two pseudoscalars has new 
contributions in the NMSSM, making it easier to obtain a sufficiently large annihilation cross sec¬ 
tion for the Higgs-pseudoscalar channel. The NMSSM has two further advantages. Firstly, we will 
consider a superpotential with only dimensionless couplings. The //-term is then dynamically gen¬ 
erated when the singlet scalar obtains a vacuum expectation value (VEV), solving the //-problem. 
Secondly, the superpotential coupling between the Higgs and the singlet gives an additional con¬ 
tribution to the Higgs quartic coupling that raises the Higgs mass at tree level. This alleviates the 
need for large stop-sector soft masses to increase the Higgs mass, improving the naturalness of the 
model. In order to obtain the observed Higgs mass, the Higgs-singlet coupling needs to be quite 
large, A > 0.7. We will focus on this regime, although smaller values of A can also provide different 
accounts of the GCE (at the expense of naturalness). 

Generically, the lightest neutralino in the NMSSM is an admixture of electroweak gauginos, 
Higgsinos and the singlino. Higgsino- and wino-dominated lightest supersymmetric particles (LSPs) 
annihilate too efficiently to be produced as thermal relics via freeze-out in the mass range of interest 
for the GCE, leaving only bino- and singlino-dominated compositions as viable. The Higgsino 
fractions of these species are strongly constrained by direct-detection experiments. For light enough 
neutralinos, these are similarly constrained by the invisible decay width of the Higgs. We will derive 
approximate analytical expressions for the resulting bounds on the Higgsino fractions. In order to 
explore the region of parameter space where the NMSSM can explain the GCE, we have performed 
numerical scans using the programs NMSSMTools 4.4.0 [32-34] and micrOMEGAs 3.0 [34]. For the 
fofo-channel, we focus on bino-dominated neutralinos. Anticipating that an enhancement of the cross 
section is required to generate a sufficiently large 7 -ray flux, the search is optimised for pseudoscalar 

1 The ink-channel was recently studied in the context of the MSSM in Ref. [26]. 
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mediators which are close to resonance in the s-channel. For the Higgs-pseudoscalar channel, we 
consider both bino- and singlino-dominated neutralinos. Our scans find regions of parameter space 
where the GCE can be generated via these channels. We study how these regions will be tested 
by upcoming direct-detection experiments and improvements in the measurement of the invisible 
Higgs decay width. 

We also consider a model that gives rise to DM annihilation into pairs of CP-even or CP-odd 
scalars. A simple extension of the NMSSM is the introduction of right-handed neutrino superfields 
and a scale-invariant term in the superpotential that couples them to the singlet [35]. We first 
consider the lightest right-handed sneutrino as the DM candidate (the corresponding cross sections 
are no longer p -wave suppressed because the DM is now a scalar particle). The superpotential term 
gives rise to a four-point scalar coupling that allows for the direct annihilation of the sneutrino into 
two Higgs or Higgs-sector pseudoscalars, in addition to s- and t/u -channel mediated processes. We 
find that the annihilation into two Higgs is typically subdominant compared to two pseudoscalars. 
We have performed a fit to the GCE for the latter final state, and find that the fit is again improved 
for pseudoscalar masses below the Higgs mass compared to the annihilation channel to Higgs-pairs. 
The sneutrino also offers the possibility of annihilation into light quarks as has been discussed in 
Ref. [36] but we will not consider this case further. Our work complements the recent analysis in 
Ref. [37], which performed a scan in a similar right-handed neutrino model, but for smaller values 
of A. 

Furthermore, we consider the possibility that the neutralino is the LSP and the sneutrino is 
the next-to-LSP (NLSP). As mentioned previously, the neutralino annihilation needs to occur on 
resonance for the 66 -channel in order to obtain a sufficiently large annihilation cross section. Due 
to thermal broadening of the resonance at high temperatures, the cross sections at freeze-out and 
in our galaxy are often vastly different. For many points in parameter space, the abundance of the 
neutralino is then too small to account for the DM. Due to its larger mass and interactions, the 
sneutrino freezes out before the neutralino. If the sneutrino is metastable and decays long after 
the neutralino freezes out, then the sneutrino decays can replenish the neutralino abundance. This 
decouples the generation of the observed DM abundance from the explanation of the GCE. Points 
in parameter space with an otherwise underabundant LSP may therefore become viable once these 
NLSP decays are included in the model. 

The WIMP annihilation in our non-minimal supersymmetric model includes the new Higgs- 
pseudoscalar and pseudoscalar-pseudoscalar channels that can provide just as good a fit to the 
GCE as the Higgs- and 66 -channels. Both the Higgs and the pseudoscalar can decay into two 
photons. Since they need to be produced close to threshold in order to give a good fit, this may 
lead to detectable 7 -ray lines in the spectrum (as previously described in Refs. [11; 12; 37]). The 
branching fraction of the pseudoscalar into two photons is suppressed compared to that for the 
Higgs. Although we find that the 7 -ray line from pseudoscalar decays is distinguishable from the 
continuum only in an optimal-case scenario, more sensitive 7 -ray experiments in the future may 
be able to detect it. Searching for these lines in the 7 -ray spectrum could provide a smoking-gun 
signal of the new channels. 

The paper is organised as follows. In Sec. 2 we summarise the different analyses of the GCE 
and perform our own fit for the Higgs-pseudoscalar and pseudoscalar-pseudoscalar channels. After 
a discussion of the LSP composition and suitable mediators for the DM annihilation in Sec. 3, we 
present the results of our parameter scans in Sec. 4. In Sec. 5 we introduce an extension of the 
NMSSM with right-handed neutrinos. We first discuss how sneutrino LSPs can explain the GCE 
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and then consider an interesting non-thermal production mechanism of neutralino LSPs from late 
decays of sneutrino NLSPs. Finally, we conclude in Sec. 6 and list the expressions of a few relevant 
couplings in the Appendix. 

2 The 7 -ray excess from the galactic centre 

Several analyses [1-10] of the Fermi-LAT data have found an excess in 7 -rays over known astro- 
physical backgrounds which originates from the galactic centre and peaks at energies 1-3 GeV. This 
excess has recently also been confirmed by the Fermi-LAT collaboration [13]. 

There have been several hypothesised explanations for this excess which have a purely astro- 
physical origin. One possible source is a previously unidentified population of millisecond pulsars in 
the galactic centre. The resulting spectrum is described by a power-law with an exponential cutoff 
and was found to fit the GCE well in Refs. [7; 8 ; 38]. However, this explanation was disfavoured in 
the more recent analysis of Ref. [10] compared to DM annihilation (see below) or a source with a 
broken power-law spectrum. Furthermore, Ref. [39] found that it is difficult to explain the spatial 
extent of the GCE with a population of millisecond pulsars while Ref. [40] argued that they can 
only account for 5% of the GCE. Another proposed explanation is that the signal originates from 
a transient cosmic-ray injection [41; 42]. Even though astrophysics may ultimately be responsible 
for the excess, we will not consider these possibilities further. 

We will instead consider the exciting possibility that the GCE is produced by the annihilation 
of DM. As we have mentioned in the introduction, we will focus on the annihilation channels of 
66 -pairs, Higgs-pseudoscalar pairs and pseudoscalar pairs, where the pseudoscalars subsequently 
decay into 66 -pairs. The latter two channels arise quite naturally in the NMSSM and the NMSSM 
with right-handed neutrinos, respectively. In Table 1, we summarise the best-fit values for the DM 
mass and annihilation cross section for the 66 -channel from different analyses 2 of the Fermi-LAT 
data (see [1-4] for earlier studies). Here and below, we denote the annihilation cross section to 
particles x and y in our galaxy (i.e. at small velocities) by ( cr xy v)o■ Note that assumptions about 
the parameters of the DM halo profile differ for the different analyses and thus the masses and 
cross sections cannot be directly compared. These values are nevertheless indicative of the target 
masses and cross sections that DM models must satisfy in order to explain the GCE. 

We next discuss the annihilation channels to a Higgs-pseudoscalar pair and two pseudoscalars 
in the context of a two-Higgs-doublet model of type II with an arbitrary number of additional 
singlets. This includes the NMSSM which we focus on later. Annihilations into scalar-pseudoscalar 
pairs from the Higgs sector of the generalised NMSSM were previously considered in Ref. [23], 
where the scalar was singlet-dominated. We will instead be interested in the case where the scalar 
is the Higgs which we will denote as h. We fix its mass to 125 GeV and its couplings to SM values, 
as implemented in PYTHIA 8.201 [46]. We denote the pseudoscalar as a and set tan/3 = 3 to fix 
its couplings (which is within the region of small tan/3 that we focus on later). Furthermore, we 
assume that neither the a nor the h can decay into other scalars. We have then performed our own 
fits to the reduced spectra of [10]. For completeness, we have also performed fits for the 66 - and 

2 The inner-galaxy analysis in the original version of Ref. [9] was subject to an error. In the updated analysis, the 
best-fit mass and cross-section for the 66 -channel shift by ~ 10% to 40 GeV and 1.8 x 10 _ 2 o cm 3 /s, respectively, with 
7 = 1.28 [43]. Since the shift is small and the updated best-fit regions are still preliminary, we will list the regions 
from the original version of [9] in Table 1. 
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Table 1: Results from different fits to the GCE assuming a 100% branching fraction to 66-pairs 
and the generalised NFW profile [44; 45] (see Eq. (2.1)) with the listed parameters. We only give 
results for the “inner galaxy” analysis of [9]. Caution should be taken in comparing results as 
different methods have been employed in each study for background modelling, data selection and 
the combining of uncertainties. Note that 7 is fit to the spectrum in [8], [9] and [10], while it has 
been simply fixed in the earlier studies. For [10], however, 7 listed here is the fixed value for which 
the fit to the mass and cross-section has been performed, while the authors find that 7 = 1.28 
provides a better fit for fixed m DM = 49 GeV. Uncertainties listed here for [6], [7] and [10] do not 
include those arising from the DM halo distribution. 


hh- channels. We assume a generalised NFW profile [44; 45] 
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with slope parameter 7 = 1.26, scale radius R s = 20 kpc and the DM density p© = 0.4GeV/cm 3 
at the radial distance of the sun from the galactic centre r©. We use the prompt photon spectrum, 
dNy/dE, for annihilations into bb from PPPC4MID [47] including electroweak corrections [48] and find 
good agreement with our own simulation using PYTHIA 8.201 [46]. We only consider the dominant 
decay channels of the pseudoscalar, to 66, t + t~ , cc, photons and gluons, and simulate the resulting 
prompt photon spectra for hh, ha , and aa final states using PYTHIA 8.201 [46]. Note that these 
spectra are unaffected by a possible singlet admixture of the Higgs or pseudoscalar. Indeed, such 
an admixture reduces their total decay widths but the branching fractions remain unaffected (to 
leading order). The differential flux measured by Fermi-LAT is given by 


dN 

~dE 


(crv)o dN^j 
8 ™dm dE 



( 2 . 2 ) 


with the line-of-sight integral fi os ds over the squared DM density. The coordinate r is centred 
on the galactic centre and can be expressed as r 2 (s,^) = r© + s 2 — 2r©scos^, where s is the 
line-of-sight distance and i/j is the aperture angle between the axis connecting the earth with the 
galactic centre and the line-of-sight. If DM annihilates into multiple final states, the different fluxes 
are summed over. 

We use the reduced spectrum of the GCE from Ref. [10] and the corresponding covariance matrix 
of the flux uncertainties including statistical and systematic errors which is publicly available. For 
simplicity we keep the pseudoscalar mass fixed and perform a two-parameter fit in the DM mass 
and annihilation cross section. In Figs, la and lc, we show the resulting best-fit spectra from 
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DM annihilation together with the spectrum of the GCE for different values of the pseudoscalar 
mass and the different annihilation channels. The smallest pseudoscalar mass in these figures is 
chosen such that Higgs decays to pseudoscalars are kinematically forbidden, i.e. m a > mft/2. The 
salmon-colored boxes depict the empirical model systematics as described in Sec. 4.2 of Ref. [10], 
the error bars correspond to the statistical errors, and the yellow boxes are the combination of 
the statistical errors, empirical model systematics and other systematics modelled as dN/dE res = 
6 x 10 -8 GeV _ 1 cm _ 2 s _ 1 sr _:L (E'/lGeV) _3 . See Sec. 4 in Ref. [10] for more details. 

Notice that, similar to the hh- channel, the spectra for the ha-channel in Fig. la have a peak 
at energies rrih /2 which is produced from on-shell decays of the Higgs to two photons. The peak is 
less pronounced than for the hh -channel because there is only one Higgs in the final state. There 
is only mild line broadening for the best-fit masses, since the hh- or ha-pair is produced close to 
threshold (see Fig. lb discussed below). If the DM mass is somewhat above threshold, however, 
the line broadening efficiently smears out the peak. Notice that the spectra have no visible peak at 
energies m a /2 from pseudoscalar decays to two photons. This can be understood as follows: In a 
two-Higgs-doublet model, the coupling of top quarks to the pseudoscalar compared to a Higgs with 
SM couplings is suppressed for tan/3 > 1 as g a tt/9h SM tt = 1/tan/3, whereas for bottom quarks it is 
enhanced as g a bb/9h SM bb = ta n /3- For intermediate values of tan/3, tan/3 < 5, the top-quark loop 
dominates the pseudoscalar decay to two photons. For the Higgs, neglecting interference effects for 
simplicity, the relative strength of the top-quark loop compared to the dominant W-boson loop is 
~ 1/10. Given that decays into bottom-quarks dominate the pseudoscalar decay width, we can then 
estimate the ratio of branching fractions for the decays into two photons of the pseudoscalar and the 
Higgs as Br(a -> 77 )/Br(h SM -r 77 ) 1/(10 x tan 4 /3). Even modest values of tan/3 (we consider 

tan/3 = 3) result in the peak being undetectable. Indeed, the peak in Fig. la for the ha-channel 
is purely due to Higgs decays. Using the Fermi-LAT limits on 7 -ray lines [49], Ref. [10] estimated 
that the line strength from the decay hh —> 4y is just below current limits and may be detected 
(or excluded) in the near future. In contrast, the intensity of the line produced by pseudoscalar 
decays is too weak to be in tension with line searches. We find that only in the best-case scenario 
of m a ~ 150 GeV and tan (5 ~ 1 is the peak barely distinguishable from the continuum. More 
sensitive 7 -ray experiments may be able to detect the photon peak for lower m a and larger tan /3 
in the future. The line at rrih /2 from Higgs decays may allow different channels to be distinguished 
because it has a different height relative to the continuum for the annihilation channels to Higgs and 
Higgs-pseudoscalar pairs and vanishes altogether for the pseudoscalar-pseudoscalar and 66 -channels. 

In Figs, lb and Id, we show the best-fit regions in the DM mass and cross section for the ha- 
and aa-channel, respectively. For comparison, we also show the best-fit regions that we find for 
the 66 - and 6 , 6 -channels. These regions agree well with those in Ref. [10]. Similar to Refs. [11; 12], 
we will assume a multiplicative astrophysical-uncertainty factor A for our best-fit cross-sections. 
This factor takes into account the uncertainties in the local DM density p©, the scale radius R s 
and the slope parameter 7 . We have used the same reference values for p©, R s and 7 for our fits as 
Ref. [11] and will therefore use their estimate for the range of the astrophysical-uncertainty factor, 
A 6 [0.17,5.3]. Notice from Figs, lb and Id that the best-fit regions for the ha- and aa-channel 
lie very close to threshold (except for the case m a = 63GeV for the aa-channel). The best-fit 
cross sections are fairly independent of the pseudoscalar and DM mass for the 6 ,a-channel, while 
the cross sections tend to increase with the pseudoscalar mass for the aa-channel. In Table 2, we 
show the values and la-regions for the best-fit DM mass and cross section and the associated y 2 
and p-values for the different annihilation channels. Notice that smaller pseudoscalar masses lead 
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channel 
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63 
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ha 

90 
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q 0+O.8 
—0.7 

25.3 

30.3 

0.28 
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137.51^0 

4 n+0-8 
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36.0 

0.03 


63 

73.0+}^ 

2 6 +0,5 
z -°-0.5 

24.3 

0.33 

aa 

90 

120 

90.0^° 0 9 

120 - 01 ^ 

q 1 +0.4 
°* 1 -0.7 

0 c+0.7 

U.O_o.7 

24.4 

31.0 

0.33 

0.10 


150 

150.01^ 

0 q+0.8 
°- y -0.8 

41.4 

0.01 


Table 2: Best-fit values and la-regions from our fits for the bb, hh, ha and aa annihilation channels, 
assuming a SM-like Higgs and tan B = 3 for pseudoscalars in the final state. 


to a better fit, i.e. a smaller y 2 and thus a larger p-value. In particular, for m a < 120 GeV, the fit 
for the ha- and aa-channels is better than for hh final states and for m a = 63 GeV (slightly above 
rrih/ 2), the fit becomes better than for bb final states. Furthermore, notice that the best-fit cross 
sections for the scalar channels are larger than for the 66 -channel. This is because the heavier DM 
requires larger cross sections to obtain the same photon flux due to the smaller number density 
(cf. Eq. (2.2); this is partly counteracted by the larger number of photons per decay since the scalar 
channels dominantly produce four 6 -quarks). 

Constraints on the DM annihilation cross section can be derived from 7 -ray observations of 
dwarf spheroidal galaxies, which are relatively clean of astrophysical backgrounds compared to the 
Milky Way. For DM which annihilates with a 100% branching fraction into 66 -pairs, the Ferrni- 
LAT collaboration finds the limits (a^v )0 < 5.00(7.91) x 10~ 26 cm 3 /s for m DM = 25(50) GeV 
[50] (with the exclusion limit growing with increasing DM masses; see their Table VI). The cross 
sections in Table 1 are consistent with these upper bounds. Recently, the Fermi-LAT collaboration 
has presented updated limits from dwarf spheroidals that would rule out cross sections {a b yv )0 > 
1.5 x 10 ~ 26 cm 3 /s for m DM ~ 30 — 50 GeV [51]. However, given that these limits are still preliminary, 
we will impose the limits from [50] but will also illustrate how our findings are affected if the updated 
limits are imposed as well. 

There is no dedicated analysis for annihilation into hh, ha and aa for either the published nor 
the latest, preliminary dwarf-spheroidal searches [50; 51]. However, 6 -quarks are still the dominant 
decay product for these channels. We can therefore use the limits for the 66 -channel by observing 
that DM of mass ?n DM annihilating into four 6 -quarks behaves kinematically similar to DM of 
mass m DM /2 annihilating into two 6 -quarks. We assume a 100% branching fraction for the decay 
h, a —> bb. Using the preliminary limit for the 66 -channel for a 50 GeV DM mass from [51] and 
accounting for a factor of 2 due to the smaller DM density and larger number of 66 -pairs in our case, 
we obtain the rough estimate {av )0 < 3 x 10~ 26 cm 3 /s for DM with m DM ~ 100 GeV annihilating 
into hh or aa close to threshold. Using a more sophisticated argument, we have checked that this 
also applies for the ha channel and depends relatively weakly on the pseudoscalar mass in the range 
of interest here. For heavier DM masses the bound on {av )0 is weaker. 
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(a) 


(b) 




(c) (d) 

Figure 1: (a) Spectrum of the GCE from Ref. [10] and best-fit spectra for DM annihilation to a 
pseudoscalar and a Higgs for different values of the pseudoscalar mass. For comparison, we also 
show the spectra for the bb- and h/i-charmels. The salmon-colored boxes depict the empirical model 
systematics as described in Sec. 4.2 of Ref. [10], the error bars correspond to the statistical errors, 
and the yellow boxes are the combination of the statistical errors, empirical model systematics and 
other systematics modelled as dN/dE ies = 6 x 10~ 8 GeV _1 cm _2 s _1 sr _1 (E/lGeV) _3 . See Sec. 4 
in Ref. [10] for a detailed discussion, (b) Best-fit regions in the DM mass and cross section for 
m a = 63,90,120 and 150 GeV (from left to right). For comparison, we also show the best-fit 
regions for the bb- and hh -channel with rrih = 125 GeV. In all regions, the inner, middle and outer 
contours delimit the lcr-, 2a- and 3cr-regions, respectively. Figures (c) and (d) are the same as (a) 
and (b) but for DM annihilation to two pseudoscalars. 
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The DM interpretation of the GCE is also constrained by other astrophysical observations, in 
particular by measurements of the anti-proton flux and radio signals from the galactic centre. Ref¬ 
erence [52] found that the corresponding limits are in strong tension with the required annihilation 
cross section for the cc- and 66-channel. However, it was subsequently argued that systematic un¬ 
certainties in the modelling of both the expected anti-proton flux [53] and the radio signals [54] were 
underestimated in Ref. [52]. Taking these uncertainties into account, the required cross sections 
may still be consistent. 

For the 66-channel, we will study DM particles with masses in the range m DM G [30, 50] GeV 
and annihilation cross sections to 66-pairs in our galaxy in the range (<7^)0 G [0.5, 5] x 10~ 26 cm 3 /s. 
These ranges are motivated by the best-fit regions in Tables 1 and 2 after taking the astrophysical- 
uncertainty factor A into account. The upper limit on the cross section is determined by the 
dwarf spheroidal limit [50]. For the 6,a-channel, on the other hand, we will consider masses m DM G 
[jm/,, 150 GeV] and require that m a + m,h < 2 m DM < m a + m h + 40 GeV, as well as cross 
sections ( cr a hv)o G [0.4,30] x 10 _26 cm 3 /s. Finally, for the aa-channel, we allow for masses m DM G 
[\rri} j, 150 GeV] with m a < ?n DM < m a -|-30 GeV and cross sections (cr aa v) o G [0.3, 30] x 1CV 26 cm 3 /s. 
These regions are obtained from Table 2, again taking the possible range of A into account. We 
will comment below on how the upper limits of the cross sections for the ha- and aa-channels are 
affected by our estimate of the dwarf spheroidal limit (the vast majority of points from our scan 
will be found to lie below this). 

3 Implications of the 7 -ray excess for SUSY 
3.1 The natural NMSSM 

The absence of new particles beyond the SM in collider searches is undermining the naturalness 
of the MSSM. The recent discovery [55; 56] of the Higgs with mass mh ~ 125 GeV and SM- 
like couplings poses additional challenges for naturalness. 3 Indeed, in order to reproduce the 
measured value for the Higgs mass in the MSSM, large soft-breaking masses in the stop sector, 
|mQ 3 |, \m Ui \, \A t \ 3> v, are necessary. These raise the quartic Higgs coupling in the effective Higgs 
potential and thus increase the Higgs mass. However, the same soft-breaking masses provide large 
corrections to the Higgs-sector mass terms via renormalization group running. The soft parame¬ 
ters at the messenger scale then have to be fine-tuned in order to stabilize the Higgs VEV at the 
electroweak scale. This tuning is reaching a level of 1 in 1000 or worse in the MSSM (see e.g. [62]). 
Much of the viable parameter space of the MSSM thus faces the little hierarchy problem. 

Naturalness is improved in the NMSSM [63-65] which adds a singlet superfield S to the particle 
content of the MSSM. We will assume the superpotential 

W D A SH u H d + ^S 3 (3.1) 

for the Higgs sector, where H u and Hd are the two Higgs doublets of the MSSM. The coupling A 
can be chosen real and positive by appropriate field redefinitions, whereas n is real but can have 
both signs. Indeed, the superpotential term A SH u H<i gives an additional contribution of order A 2 to 
the Higgs quartic coupling. This raises the Higgs mass at tree-level, circumventing the necessity of 
heavy stop-sector soft masses and thus relieving the tuning. In addition, the A 2 -contribution to the 

3 See [29; 57-61] for the connection between the Higgs couplings and naturalness. 
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quartic coupling makes the Higgs VEV less sensitive to destabilising loop corrections, as follows from 
the minimisation conditions for the Higgs potential [66]. Due to this effect, the required amount 
of tuning for the Higgs VEV decreases with increasing A. It was pointed out in [67], however, that 
the Higgs mass in turn grows and eventually overshoots the observed value for increasing A. Loop 
corrections and mixing effects are then required to bring the Higgs mass back down. This is an 
additional type of tuning that should be taken into account when assessing the naturalness of the 
model and a combined tuning measure was defined in [67] as the product of the tunings in the 
Higgs VEV and the Higgs mass. This combined tuning was found to be minimised for intermediate 
values A r\_/ ~\_ m 

The superpotential (3.1) does not have any dimensionful parameters. The //-term is instead 
dynamically generated from the term A SH u Hd once the scalar component of S obtains a VEV. 
This provides a solution to the //-problem of the MSSM. A parameter scan of this model was 
performed in Ref. [67], with a focus on regions with A ~ 1 in order to minimise the fine-tuning. 
Two model-building assumptions were made to improve naturalness further. The first of these is 
a split spectrum, where the first-and-second generation squarks and all sleptons are much heavier 
than the remaining superpartners, which allows collider constraints to be satisfied while keeping 
the superpartners most relevant for naturalness relatively light. Secondly, a small messenger scale 
of 0(20 TeV) reduces the size of loop corrections during the renormalization group evolution and 
therefore the required amount of fine-tuning. Such a low scale is further motivated by the fact 
that, for A > 0.7, the coupling runs fast enough to reach a Landau pole below the GUT scale. 
Some UV completion for the Higgs sector is then required to emerge below this energy scale. In 
order to ensure that the low-energy description of the theory is valid up to the messenger scale, 
one has to demand that no Landau pole appears at lower energies. Since the energy scale at which 
the Landau pole arises decreases with increasing A, a low messenger scale allows for a larger range 
of A. Imposing the relevant limits on superpartner masses, electroweak precision tests and flavour 
violation, Ref. [67] found large regions of parameter space for which the combined tuning in the 
Higgs VEV and Higgs mass is still better than 5%. These regions had parameters in the ranges 
0-8 < A < 1.3, 2.5 < tan/3 < 4.2, and 100GeV < |// e ff| < 600GeV, where // e fr = Xv s is the effective 
//-term after the scalar component of S obtains a VEV v s . 

Motivated by naturalness, we focus on the region A ~ 1 for the NMSSM. In order to avoid the 
//-problem, we consider the scale-invariant superpotential (3.1). For definiteness, we furthermore 
assume a split sparticle spectrum. We discuss this spectrum in more detail in Sec. 4. We emphasise, 
however, that our results generically hold for non-split spectra. The messenger scale, on the other 
hand, does not play a role in our analysis because we specify all parameters near the electroweak 
scale. 


3.2 Composition of the LSP 


We will focus on neutralinos as the LSP. The symmetric neutralino mass matrix in the basis 
(B,W 3 ,h° d ,h° u ,~s) is given by 


/ M x 


M * 


0 

M 2 


— cos (3 

^v sin/3 

0 

-%ucos/3 

\j i 

—-%usin/3 

\j i 

0 

0 

/^eff 

—Xv sin j3 


0 

—Xv cos (3 



2kv s 


(3.2) 
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where M\, M 2 are the gaugino masses and gi,g 2 are the gauge couplings corresponding to U(l)y and 
SU(2)i, respectively. We denote the neutral scalar components of the superfields H u , Hj and S by 

hi, h° d and s, respectively. The electroweak scale is then given by v = \J (hi ) 2 + (h d ) 2 = 174 GeV, 
while tan/3 = (hi) / (h d ) and v s = ( s ). 

A Higgsino- or wino-dominated neutralino with the required properties for the DM and to 
account for the GCE is problematic. The LEP bounds on charginos require that M 2 , g e s > 
100 GeV [ 68 ], which restricts the possibility of a Higgsino- or wino-dominated LSP in the mass 
range relevant for explaining the GCE via the 66 -channel. Furthermore, neutralinos with a signif¬ 
icant admixture of Higgsino states are severely constrained by direct detection [69; 70]. Although 
very pure Higgsino-like neutralinos avoid these limits, in our mass range of interest they are in 
tension with indirect-detection bounds [70] and have annihilation cross-sections that are too large 
for thermal production (the latter being generically true for Higgsinos of mass < ITeV). This 
problem also arises for winos with masses below 2.7 TeV, which are also in strong tension with 
indirect-detection bounds [70-72]. This leaves bino- or singlino-dominated LSPs. For simplicity, 
we will decouple the wino, as it need not be lighter than several TeV in a natural superpartner 
spectrum. 

Both a bino- and a singlino-dominated LSP will generically have a non-negligible Higgsino 
admixture, while the admixture between the bino and the singlino is naturally suppressed, as 
follows from the mass matrix (3.2). Both direct-detection experiments and constraints on the 
invisible decay width of the Higgs restrict the size of the Higgsino admixture. In the scale-invariant 
NMSSM, these bounds are readily satisfied if g e g = Xv s is sufficiently large to suppress the Higgsino 
mixing. For the bino-dominated case, the mixing matrix elements with Higgsinos in (3.2) are 
relatively suppressed by the gauge coupling gi, so smaller values of g e g are more tolerable. For the 
singlino-dominated case, however, the Higgsino mixing is greater because A ~ 1. The latter case 
therefore prefers a larger g, e s, which in turn increases the fine-tuning if too large. Furthermore, 
simultaneously keeping the singlino mass (the (5, 5)-element in (3.2)) small, and hence the singlino- 
dominated LSP light, requires a relatively small singlet self-coupling k ~ 0.1. These features make 
the parameter space with a viable singlino-like LSP more restrictive. For the bino-dominated LSP, 
the direct-detection bounds (in conjunction with the other experimental constraints on the model) 
are most easily accommodated in the region of g e Q < 0, as observed previously in [70]. 

To obtain the constraints on the Higgsino fraction from direct detection and the invisible decay 
width of the Higgs, we first derive the resulting bounds on the coupling c/,^ of the Higgs to LSPs. 
The coupling is given by [73] 


°hxx ~ 5 i Mi (S11A/13 - S12A/14) — 52 A/12 ( 5 nA/l 3 — 5 i 2 A/14) 

— y/2 k S 13 A/”ig + \/2A (5n J\T\a A /15 + £12 A /13 A /15 + £13 A /13 Mu) , (3.3) 


where S and M are diagonalisation matrices that relate the mass eigenstates hi and x.i of the 
CP-even Higgses and the neutralinos, respectively, to the interaction eigenstates according to 



(x A 


f B \ 

X2 


w 3 

X3 

= M 

~K 

X4 



\X5/ 


\ s ) 


(3.4) 
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We identify the lightest CP-even state h\ with the Higgs observed at the LHC. For simplicity, this 
state will be abbreviated as h and the LSP xi as X- 

The Higgs contributes to spin-independent DM-nucleon scattering. Electroweak gauge bosons 
and pseudoscalars, on the other hand, contribute only to spin-dependent scattering for which 
the limits are much less stringent. The direct-detection experiment LUX provides the strongest 
constraints on spin-independent DM-nucleon scattering to date (which are also strongest in the mass 
range relevant for the GCE) [74], For DM masses ~ 40 (200) GeV, the bound is u® 1 < 7 (25) x 10~ 46 
cm 2 . In the absence of light scalars that could interfere with LSP-nucleon scattering (coloured 
sparticles or the other CP-even states from the Higgs sector), this directly provides a bound on the 
coupling Chxx■ The spin-independent cross section for LSP-nucleon scattering is given by (assuming 
isospin symmetry for simplicity) [75] 


(T 


SI 


^ 2 r 2 

- - /Lv/p , 

7r 


where fijsi ~ m p is the LSP-nucleon reduced mass and 


(3.5) 


f p = m p (^2 f: 


q=u,d,s 


f 'P °XL 

Tq m n 


4- — f p 

+ 27 j t g 


X 

q=c,b,t 


m n 


(3.6) 


is the proton form-factor with f p = {0.0153,0.0191,0.0447} and /{^ = 1 — /£ — — fj, , while 

m p is the proton mass and m q denotes the quark masses. For tree-level parton scattering mediated 
only by the Higgs h, 


m u Sl2 Chxx j _ m d *^11 c h\X 

2\[2v sin /3 ml d 2\/2ucos /3 ml 


(3.7) 


and analogously for the other up- and down-type quarks. For a Higgs with SM-like couplings 
(as observed at the LHC), we have <Sn ~ cos (3 and 5i2 ~ sin/3. The bound on the cross-section 
from LUX then gives the approximate bound | c *-xxI ~ 0.04(0.07) for m^ = 40 (200) GeV. In this 
derivation, we are crucially assuming that the other CP-even Higgses are sufficiently heavy to 
neglect their contribution to the scattering cross-section. This is justified for most of the points in 
our analysis and was found to hold over much of the space uncovered by Ref. [67]. 

If Higgs decays to LSPs are kinematically allowed, the coupling c/^ is also constrained by 
limits on the invisible branching fraction of the Higgs provided by fits of the Higgs couplings to 
LHC data. These require that Br(/i —> inv) < 0.24 [76; 77], taking the fit to the Higgs couplings 
with assumptions of SM-couplings to gluons and photons and couplings to weak vector bosons of 
at most SM-strength (the latter assumption being always satisfied by a general two-Higgs-doublet 
model with additional singlets as in the case of the NMSSM). The decay width of the Higgs to 
LSPs is given by 


r h 


—>xx 


\°hxx\ 


2 m h 
16-7T 



(3.8) 


at tree-level. For a Higgs with SM-like couplings, the total decay width to SM particles is Tsm ~ 
4.07MeV [78]. For m ^ = 40 GeV, the bound on the invisible branching fraction then gives | c fexxI ~ 
0.03 which is slightly more stringent than the bound from direct detection, although this weakens for 
masses closer to threshold (for = 50 GeV, \chxx\ 0.04). As DM-nucleon scattering mediated 
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by the Higgs may be subject to interference from amplitudes mediated by other light scalars, the 
bound from the invisible branching fraction is more robust provided that the DM is light enough. 

These limits on the coupling c^x translate into limits on the Higgsino fraction of the LSP 
under simplifying assumptions. Let us define c max = 0.03 (0.07) as the maximum allowed Higgs- 
LSP coupling for m x = 40 (200) GeV. We first consider a bino-dominated LSP. From the neutralino 
mass matrix (3.2), we see that the mixing matrix element between the bino and singlino vanishes 
and thus A /15 “C 1, while we assume that the wino is decoupled so that A /12 1. The Higgs-LSP 

coupling is then approximately given by 

c hxx ~ <7i A/11 (S11 A/13 — «Si2 A/14) + \/ 2 A S13 A/13 A/14 ■ ( 3 . 9 ) 

Anticipating a small Higgsino fraction, we set A /11 ~ 1- Furthermore, we use 5n ~ cos /3 and 
5,2 ~ sin/3 for a SM-like Higgs. Then ignoring the possibility of blind spots and applying the 
limits on || to the first and second term in (3.9) separately gives the approximate bounds 
IA/ 13 I < C max /(51 cos/3) and |A7 i 4 | < C max / (51 sin/3). 4 

Let us now consider a singlino-dominated LSP. Since the mass matrix elements mixing the 
gauginos and the singlino vanish, we have A/ 11 , A /12 1 . The Higgs-LSP coupling then simplifies 

to 

c hxx ~ V^A ( 5 n A/14 A/15 + 5 i 2 A/13 A/15 + S13A/13A/14) — a/2k5i3 A/^5 . ( 3 . 10 ) 

Anticipating a small Higgsino fraction, we set A /15 ~ 1- We again ignore possible blind spots and 
apply the limits on c^x to the first two terms in (3.10) separately. This yields the approximate 
bounds |A/i 31 < c max /(\/2Asin/3) and |A/i 4 | < c max /(\ / 2 Acos /3 ). 5 

We caution the reader that the above analysis is quite simplified and that the bounds derived 
are only approximately valid. In particular, blind spots, as defined in [70], can potentially be 
important and open regions of the natural parameter space that would otherwise be ruled-out. 
These involve the partial cancellation of the terms in the couplings and allow for larger Higgsino 
components of the neutralino than naively derived above. Furthermore, a singlino-dominated LSP 
requires small k in the scale-invariant NMSSM and in this region a CP-even singlet-dominated 
state from the Higgs sector is typically light and can contribute to LSP-nucleon scattering, thereby 
complicating the analysis. However, we still find that our estimates above are approximately 
consistent with the results from our scans below. For the singlino-dominated LSP, we find that 
typically |A/ 14 1 ~ 0.3 — 0.4 and |A/ 13 1 ~ 0.05 for ~ 100 GeV, while for the bino-dominated LSP, 
IA/ 14 1 < 0.01 and |A/i 3 1 ~ 0.1 ( 0 . 2 ) for ~ 50 (100) GeV. Given the range of A and tan/3 that we 
scan over, these values are approximately consistent with our bounds (although A /13 can be up to 
a factor of ~ 2 larger for the singlinos). 

4 Using these bounds, the third term in (3.9) satisfies e.g. \/2\\ S 13 A /13 A/ 14 1 < 0.7 c ma x A/ sin 2/3 for = 40 GeV. 
Since A ~ 1 and tan (3 < 4, this term at most marginally violates the limit on Chxx and can thus give only a slightly 
more stringent limit on the product A /13 A/14 than the first two terms. 

’’Ignoring possible cancellations among the terms in (3.10), the last term translates to a bound on the singlet 
fraction of the Higgs, |<Si 3 1 < c max / (\/2\k\). This is, however, not constraining in the limit of small k (as is required 
for a singlino-dominated LSP). But even for S 13 ~ 1 (which would typically be in conflict with the measured Higgs 
couplings), the third term in (3.10) satisfies e.g. \/ 2 A|Si 3 A /13 A/ 14 1 < 0.04 c max /(A sin 2/3) for = 40GeV. Since 
A ~ 1 and tan (3 < 4, this term is always negligible when deriving bounds on the Higgsino fractions. 
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3.3 Mediators for the annihilation of the LSP 


We will next consider possible mediators for the annihilation of DM. Due to the hierarchy in the 
Yukawa couplings, annihilations via the Higgs naturally give a large branching fraction to 66-pairs 
for an LSP within the best-fit mass range for the 66-channel. However, annihilations of Majorana 
fermions via scalars are p -wave suppressed and the annihilation cross section via the Higgs (and 
the other two CP-even Higgses) in our galaxy is consequently too small to account for the GCE. 
Annihilations into fermions mediated by the Z-boson are helicity-suppressed in the s-wave and are 
also too small to significantly contribute to the GCE. In the early universe, however, the velocity 
suppression of the p -wave component mediated by the Z-boson is lifted and can enhance the total 
annihilation cross section during freeze-out. 

For bino-dominated LSPs, 6-channel exchange of sbottoms is another possible annihilation pro¬ 
cess to 66-pairs. However, the cross section is typically too small because of the small hypercharge 
coupling and collider limits on the sbottom mass. It was found in Ref. [24] that a highly degenerate 
sbottom-neutralino pair with mass splitting 2 GeV < — m% < nib and large mixing in the sbot¬ 

tom sector can have a sufficiently large annihilation cross section and may avoid direct-detection 
constraints. 

As we will now discuss, CP-odd scalars from the Higgs sector are suitable mediators. The anni¬ 
hilation of Majorana fermions via these states at small velocities is dominantly s-wave and thus not 
suppressed. In addition, CP-odd scalars contribute only to spin-dependent DM-nucleon scattering 
for which the bounds are much weaker. One linear combination of the CP-odd components of the 
Higgs doublets /i° and h d gets eaten by the Z-boson. We denote the orthogonal linear combination 
of the doublets by Hj and the CP-odd component of the singlet scalar by sj. In the basis {Hj, s/), 
the mass matrix then reads 


M, 


rn. 


CP-odd 


Au(t1a — 2nv s ) 

A v 2 sin 2/5 + 2 n'j — 3 kA k v s 


(3.11) 


where m\ = 2A v s (A\ + kv s )/ sin2/5. We write the lightest CP-odd state as a = cos 6 A 4 P/+sin 6 ^ 457 . 
The coupling of a to an LSP pair then reads 


^ a xx — ^ cos 


Oa [92 X 12 (A/"i 3 sin/5 - A/j 4 cos/5) - gi A/ii(M 3 sin/5 - JVj 

-L- a/T \ \f. - ( ™ one /■? -L- \f, , cin /L -L- n a /*5 Gin Q a ( \ A f, . 




We can use the upper bounds on the Higgsino admixture of the LSP, A /13 and A/ 14 , derived in 
the last section to bound the size of the coupling c a ^y. For the bino-dominated case, using A /11 ~ 1 
and A/ 12 , A /15 < 1 as before and the triangle inequality, we find that 

\ c axx\ ~ 2 Cmax | COS 0a! /sin 2/5 + 2 V 2 A Cmax I sin @a \/(ffi sin 2/5) . (3.13) 

For the singlino-dominated case, on the other hand, again using A /15 ~ 1 and A/ 11 , A /12 1, 

|CaxxI ~ 2 c max | cos 9a |/ sin 2/5 + y/2 \ sin0A| (|«| + c^ ax /(Asin2/5)) . (3-14) 


From these bounds, we may estimate the available parameter space that can explain the GCE for 
the annihilation channels to bb and to Higgs-pseudoscalar. To this end, we will discuss these two 
channels separately. 
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Figure 2: Feynman diagrams contributing to neutralino annihilation to 66 -pairs. 


3.3.1 Annihilation into 66 


The most important processes for annihilation into 66 -pairs are shown in Fig. 2. The CP-odd 
scalars from the Higgs sector couple to the SM via the usual Yukawa couplings and can thus have 
large branching fractions to 66 -pairs. The coupling of the lightest CP-odd Higgs to a pair of bottom 
quarks is given by 


c abb ~ 



cos 9a tan /3 , 


(3.15) 


where = m^fv is the SM Yukawa coupling to bottom quarks. 

The annihilation cross section for the process XX —> a* —» bb at small velocities is given by [17] 


(av) 0 


2 x 10 


26 cm3 

c abb 

2 

C «XX 

2 / rn k \ 2 / (120 GeV ) 2 - 4(40 GeV ) 2 \ 

S 

Vb 


0.5 

V40 GeV / ( — 4?n 2 + J 


2 


(3.16) 


where r a is the total decay width of the pseudoscalar. Let us for the moment assume that the 
annihilation is not resonantly enhanced, i.e. that m a is not very close to 2 m k . The T a -dependent 
term in the denominator is then negligible and the last two factors in the formula are of order one 
in the mass range of interest for the GCE. In order to explain the GCE and avoid constraints from 
additional annihilation channels, the process XX a * bb should dominate the total annihilation 
cross section at the present time. In the absence of additional important annihilation channels 
during freeze-out (e.g. mediated by the Z-boson) and coannihilation, the cross section at freeze- 
out is of the same size as (3.16). Both the required amount of DM and the GCE can then be 
accommodated if the remaining factors in (3.16) are close to one. 

From our upper estimates of the coupling c akk , we can obtain an upper bound on these fac¬ 
tors and the cross section. Note that these estimates and the coupling c abb grow with tan /3 
(cf. Eqs. (3.13)-(3.15)) and that c abb vanishes if a consists purely of sj. In order to obtain a 
conservative bound, we therefore set tan/3 = 3, corresponding to the largest value for which we 
find points in our scan with a significant ///-admixture to a. We assume m k = 40GeV (which is 
within the mass range of interest for the 66 -channel) and distinguish between the two cases of a 
being dominantly Hj and an equal admixture of Hj and sj (the bounds become stronger if a has 
an even larger s/-admixture). Let us first consider the bino-dominated case. The second and third 
factor in (3.16) can be estimated as 


^ abb 

2 

C »XX 

2 < r o.2 

for cos 9a s 

3 1 

Vb 


0.5 

~ [(0.2 + 0.07A) 2 

for cos 9a s 

s l/x /2 


(3.17) 
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Figure 3: Feynman diagrams contributing to neutralino annihilation to ha- pairs. Note that there 
is a u-channel diagram in addition to the t-channel one. 


Similarly for the singlino-dominated case, setting again tan /3 = 3, we can estimate the second and 
third factor in (3.16) as 


^abb 

2 

C “XX 

2 < fo.2 

for cos 9a s 

3 1 

Vb 


0.5 

~ \(0.2 + 3|k| + 5 x 10 _ 3 /A ) 2 

for cos 6a s 

S l/y/2 


This factor for either the bino- or singlino-dominated case needs to be not much smaller than one for 
a sufficiently large annihilation cross section (3.16). As we have discussed in the last section, k< 1 
is required in the scale-invariant NMSSM for a light singlino-dominated LSP with a small Higgsino 
fraction. This means that the ^-dependent term in (3.18) does not help. Even though the case 
cos 0,4 ~ 1 seems marginally consistent with the required cross section for the GCE when taking 
the astrophysical uncertainty factor A into account, the freeze-out cross section is then typically 
too small to generate the required amount of DM. This would require additional contributions to 
raise the freeze-out cross section, e.g. mediated by the Z-boson. In addition, we find that tan /3 < 3 
for cos 0,4 ~ 1. This is because cos 9a is correlated with tan/3. Indeed, the mass term of the singlet 
si, the (2, 2)-element of the CP-odd mass matrix (3.11), is suppressed while the mass term of the 
doublet Hi, the (1, l)-element, is enhanced by tan/3 . 6 This means that the admixture of Hi to a 
and thus cos6 a decreases with increasing tan/3. This dependence is indeed what we observe in our 
scan. For tan/3 ~ 3, we find points only up to | cos 6a\ ~ 0.75 (for a bino-like LSP; cos 6a is even 
smaller for larger tan/3). This means that the bound for the case cos 0,4 ~ 1 in (3.17) and (3.18) 
is actually more stringent. We will therefore need an enhancement of the cross section. To this 
end, we will rely on the regime m a ~ 2 m ^ in which case the cross section is resonantly enhanced 
and the last factor in (3.16) is much larger than unity. This is in contrast to Ref. [17], which 
considers smaller A and can thus tolerate large tan /3, allowing for larger couplings c ab i and c a ^ 
and a sufficiently large cross section without resonant enhancement for the bino-dominated case 
(though the smaller A is at the expense of increased fine-tuning due to the Higgs mass). 

3.3.2 Annihilation into ha 

The annihilation of neutralinos is p -wave suppressed if the final state is even under CP, such as 
for Higgs or pseudoscalar pairs [22], Therefore it is only through the Higgs-pseudoscalar channel 
that neutralino annihilation into scalars can account for the GCE. The graphs that contribute to 
the corresponding cross section are the s-channel exchange of the two pseudoscalars and the t/u- 
channel exchange of the neutralinos. These are shown in Fig. 3. As follows from Fig. lb, a good 

6 An exception occurs if the first term in the (2, 2)-element is negative and cancels against the second term in 
which case the (2, 2)-element is enhanced by tan /3 too. 
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fit to the GCE requires that the annihilation into ha happens close to threshold. Let us define the 
dimensionless ratio 5 = ( 2 m* — (m a + m/,))/2m^ so that 5 = 0 corresponds to the annihilation 
being precisely at threshold. Expanding for 6 -C 1, the cross section for the process XX ha is 
given by [ 22 ] 

Jl ( mh \ 1/2 f i _ 1/2 

8vr [mxJ V 2m xJ 

x c aah c axx _|_ 

m h (4m * - m h ) 4 

where the heaviest pseudoscalar is denoted by 02 (we neglect the small width of the a). 

From Fig. lb, we see that the best-fit regions have m * ~ 100 GeV and 5 ~ 0.1. The best-fit 
cross section, on the other hand, depends on the value of the astrophysical uncertainty factor A but 
is consistent with the freeze-out cross section ~ 2 x 10 _ 26 cm 3 /s required for the WIMP miracle. In 
the absence of resonance enhancement, additional important annihilation channels (e.g. mediated 
by the Z-boson) and coannihilation, the cross section at freeze-out is of the same size as (3.19) 
(deviations from these conditions are tolerable to some extent given the relatively large range for 
the uncertainty factor A). Assuming the values for 6, m* and the target cross section above, we 
will now show that (3.19) can be large enough to explain both the GCE and the observed DM 
abundance without relying on special enhancements such as a resonance. Considering only the 
t-channel exchange of the second neutralino X 2 , the cross-section may be approximated as 



v/5 


c a2ah c a,2XX 


- m l 2 


T 771(22 ^“2 


+ 2 


5 

£ 

k=1 


C «XXfc c hxXk 

m h + m^ k - m 


(3.19) 
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2 x 10 -26 


cm 3 

s 


( C a\'X2 c hxx 2 

V 0.2 x 0.2 


2 — 125 GeV / 5 
2 • 100 GeV - 125 GeV V 0 T 
125 GeV - 100 GeV + 150 GeV V 
125 GeV — + m^ 2 J 


100 GeV 


m< 


(3.20) 


where we set mh = 125 GeV. We have focused on the second neutralino because the LSP necessarily 
has a small Higgsino-fraction as discussed in Sec. 3.2. This limits the size of the couplings c a ^ 
and <'hxx that would appear in (3.20) for t-channel exchange of the LSP and thereby the size of 
the cross section. Since there is no limit on the Higgsino fraction of the second neutralino, the 
couplings c a xx 2 an d Chxx 2 can be much larger. Thus, if the second neutralino is relatively light, we 
see that a sufficiently large cross section to explain the GCE and the observed DM abundance can 
be obtained. 

Similarly, considering only the s-channel exchange of the pseudoscalar a, the cross section is 
given by 


(a ha v )0 « 2xl0 - 26 — 


V v / 


0.08 


2 /100 GeV \ 


m-c 



2 my, ~ 125 GeV 
2 - 100 GeV - 125 GeV 

4-100 GeV-125 GeV' 
X 1 4 m* - 125 GeV 


(3.21) 


where we set mh = 125 GeV. The coupling c aa h between the Higgs and two lightest pseudoscalars 
is dimensionful and determined by v, v s and A\ (for brevity, we do not give the expression for this 
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coupling which can be found in [73]). Since v s and A\ can be much larger than v, the suppression 
from the coupling c a x\ for a small Higgsino fraction may be overcome by the coupling c aa h. Thus, we 
again see that a sufficiently large cross section to explain the GCE and the observed DM abundance 
can be obtained. 

However, we find in our scan below that generally neither of the terms in (3.19) dominate the 
others, but that they must all be considered and can interfere with each other. When all the 
contributions are combined they allow for a sufficiently large cross section without the need for a 
resonant enhancement. 

4 Parameter scans 

In order to find parameter sets for which the NMSSM can explain the GCE via the bb- or ha- 
channel, we have performed numerical scans using the programs NMSSMTools 4.4.0 [32-34] and 
micrOMEGAs 3.0 [34], To simplify the parameter space, we decouple the first and second generation 
of squarks and all generations of sleptons, fixing their masses at 2 TeV (similar to the scan performed 
in [67]). We similarly decouple the wino, choosing 2 TeV for its mass. Due to the small Yukawa 
couplings of these sparticles and the small gauge couplings, these masses are compatible with 
naturalness. We fix the soft masses of the stops and sbottoms at 1 TeV with vanishing trilinear 
terms Af j- } and the gluino mass at 2 TeV. These values are chosen such that the collider constraints 
on their masses are fulfilled but the fine-tuning still remains approximately optimal. Note that the 
physics of the DM is relatively independent of the properties of these particles, so as long as collider 
constraints on these masses are satisfied, we do not expect our results to qualitatively change if 
these are allowed to take different values. We have performed three different scans, for (i) a bino- 
like LSP with dominant annihilation to ha, (ii) a singlino-like LSP with dominant annihilation to 
ha and (iii) a bino-like LSP with dominant annihilation to bb. 

We impose all relevant collider and flavour constraints which are implemented in NMSSMTools 
4.4.0 (without trying to explain the discrepancy from the SM in the anomalous magnetic dipole 
moment of the muon), which in particular require a Higgs at 125 GeV with SM-like couplings 
consistent with the signal strengths measured at the LHC. Most points that pass all experimental 
constraints have m a > m/,/2 and thus decays h —> aa are kinematically forbidden. Indeed, a signif¬ 
icant decay rate h —>• aa would induce large deviations from LHC measurements in the branching 
fractions of the Higgs decays into photons, 6-quarks, r-leptons, W- and Z-bosons. We only find 
a few points with m a < rrih/2 and these are all close to the threshold. Since the pseudoscalars 
can have a significant singlet-admixture, LHC searches do not significantly constrain the range of 
pseudoscalar masses that we consider [79]. For all points in the scans, the couplings A and k do 
not evolve toward a Landau pole below 20TeV. We use micrOMEGAs 3.0 to calculate the relic 
abundance of the LSP and the cross-sections relevant for direct and indirect detection. We impose 
the range 0.107 < < 0.131 when we require that the LSP abundance be consistent with the 

observed DM abundance. This range is obtained from the best-fit value for the DM abundance 
from Planck [80] after including a ±10% theory error in order to account for uncertainties in our 
calculations (the measurement error quoted by Planck is subsumed by this). 

In addition, we check constraints from electroweak precision measurements which are not im¬ 
plemented in NMSSMTools 4.4.0. These can be important for large A ~ 1 since the mixing of the 
Higgsinos with the other neutralinos then induces sizeable deviations from the custodial symmetry. 
We calculate the Peskin-Takeuchi parameters S and T using the formulas from Refs. [81-83] for 
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AA a [GeV] 

A k [GeV] 

Vez [GeV] 

Mi [GeV] 

A 

K 

tan/3 

bino-like [-50,50] 

[-100,100] 

[-300,-100] 

[60,170] 

[0.6,1.4] 

[0.1,1.6] 

[2,5] 

singlino-like [-50,50] 

[-100,100] 

[-600,-200] 

2000 

[0.6,1.4] 

[0.05,0.5] 

[2,5] 


Table 3: Parameter ranges for the two random scans for the ha-channel. The first and second line 
are for the scan optimised for bino-like and singlino-like LSPs, respectively. 


the neutralino-chargino and stop-sbottom sector and those from Refs. [84; 85] for the Higgs and 
singlet scalars. We compare the results with the 95% exclusion region determined in [86]. As was 
found previously, for A ~ 1 the neutralino/chargino-sector tends to provide the largest contribution 
to T and tan (3 3> 1 is ruled out due to custodial-symmetry breaking. 


4.1 Parameter scan for the /ia-channel 


We have performed two random scans using NMSSMTools 4.4.0 [32-34] for the Higgs-pseudoscalar 
channel, one optimised for bino-like LSPs and the other for singlino-like LSPs. The masses of 
the squarks, sleptons, gluino and wino are fixed as mentioned previously and the remaining free 
parameters are scanned over the ranges shown in Table 3. Here A\ is partly determined by the 
requirement that the singlet admixture to the Higgs is small. To discuss this in more detail, let 
us rotate from the basis into the ( h SM ,H ) basis, where only h SM obtains a VEV and thus 

couples to SM particles precisely like the SM Higgs. The admixture of the orthogonal state H and 
the singlet s to the lightest CP-even Higgs then drives its couplings away from SM values. The size 
of the admixture with s is determined to leading order by the mass term, m s of s and its mixing 
mass term, rri} ls with h SM : 


1 


r — N 2 2 i z 

£ D -m s s + m hs h SM s , 


(4.1) 


where 


m 2 = kA k v s + 4 k 2 v 2 + -^-A a u 2 sin 2/3 , 

4lV s 

rrif ls = 2X 2 v s v — sin 2/3 A v(A\ + 2kv s ) ■ 


(4.2) 

(4.3) 


The sizes of A\ and A K are limited by the requirement that the Higgs sector is non-tachyonic. 
Furthermore, as we have discussed in Sec. 3.2, small k ~ 0.1 is required for a singlino-dominated 
LSP with a small Higgsino admixture. We find in our scan that this then implies that rn s is not 
much larger than the Higgs mass. In order to suppress the singlet mixing with the Higgs and 
the resulting deviations in the Higgs couplings, we therefore fix A\ to minimise the mixing matrix 
element rri 2 s . More precisely, we set A\ = A^ + AA a , where A^ is determined by = 0 and 
A A\ is randomly chosen in the range [—50, 50] GeV. Although a bino-dominated LSP allows for 
larger values of k, we find that // e ff must be relatively small for a sufficiently large cross section 
(as suggested in (3.20), where a Higgsino-like \2 of mass m^ 2 ~ 150 GeV is preferred). This again 
results in a relatively light singlet-like state of several hundred GeV. In order to obtain a Higgs 
with SM-like couplings, we therefore also impose this condition on A\ for the bino-dominated case. 
We have checked that our findings are qualitatively unchanged after relaxing this condition. 
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Figure 4: Results from the two random scans with parameter ranges given in Table 3 optimised 
for the ha-channel. We require a branching fraction of annihilations into ha of at least 85%. The 
left panel (a) shows the annihilation cross section at the present time, {(Jh a v) o, and at freeze-out, 
( ctv)t f ■ The right panel (b) shows the spin-independent cross section for scattering of LSPs off 
nucleons versus the LSP mass. 


Results from the scans are presented in Fig. 4. All points in the scatter plots have an LSP 
abundance in the range 0.107 < D^/i 2 < 0.131, consistent with the observed DM abundance as 
discussed above. Furthermore, all points have masses my and m a and an annihilation cross section 
i a ha v ) o into ha in our galaxy within the best-fit ranges as discussed at the end of Sec. 2. We also 
demand a branching ratio to ha of at least 85% in order to have a relatively pure annihilation 
into these states. In Fig. 4a, we show a scatter plot of the annihilation cross-sections during 
freeze-out and at late times. The approximate degeneracy of these cross sections when (crhav )o ~ 
(1 — 3) x ICR 26 cm 3 /s reflects the dominance of the XX ha process and the absence of a large 
special enhancement such as a resonance or coannihilation. The late-time cross section can reach 
values of up to ~ 2 x 10~ 26 cm 3 /s for the bino-dominated LSPs, whereas it is smaller for the 
singlino-dominated LSPs. These cross sections are nevertheless consistent with the best-fit regions 
prescribed in Sec. 2 once uncertainties in the DM distribution are accounted for. Larger late¬ 
time cross sections are too big for the LSP to freeze-out with the required relic abundance. This is 
because additional processes that are p -wave suppressed in our galaxy (most notably those mediated 
by a singlet-like CP-even scalar) activate and contribute significantly in the hotter early universe. 
Furthermore, we find that the points with bino-dominated LSPs that have significantly enhanced 
cross-sections at freeze-out, notably most of the “column” clustering at {gv)t f ~ 5 x 10~ 26 cm 3 /s 
in Fig. 4a, as well as the diffuse region up to ( <jv)t f ~ 2 x 10 _25 cm 3 /s, undergo coannihilation 
with a Higgsino-like second-lightest neutralino. The larger effective cross section at freeze-out is 
partly because the aa and hh channels are not p -wave suppressed for coannihilation since the 
Pauli exclusion principle no longer applies to the initial coannihilating particles. The few points 
with (ctv)t f ~ 10~ 25 cm 3 /s have this coannihilation resonantly enhanced by a singlet-like CP- 
even scalar. Because these enhancements have a strong temperature dependence, they diminish 
as the universe cools such that the LSP still freezes-out with the observed DM relic density (the 
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cross sections plotted in Fig. 4a are calculated for the instantaneous temperature at which the 
LSP begins to freeze-out). We emphasise, however, that neither a resonance nor coannihilation is 
necessary for a bino-like LSP to explain the GCE as can be appreciated by the fact that there are 
many points with {crha v )o ~ {ov)t f in Fig. 4a. We have also plotted, as a dashed line in Fig. 4a, 
our estimate from Sec. 2 of the limit on (<Jha,v )o arising from the latest Fermi-LAT study of dwarf 
spheroidals. All our points satisfy this limit. 

The characteristics of the two types of LSPs are quite distinct. The points with a singlino-like 
LSP are dependent upon a large |A\| ~ 1 — 2.5 TeV. This is because the neutralinos x.2 and X3 
are typically relatively light and Higgsino-dominated, leading to a large contribution to the cross 
section. Although the significance of the f-charmel exchanges of X‘i and X 3 in (3.20) diminishes as 
| increases (as this determines the Higgsino masses), we only scan over relatively small values 
l/^effl < 600 GeV to try to maintain naturalness. The contribution to the cross section from the 
s-channel exchange of the a is then required to be subdominant. To achieve this, A\ is selected to 
cancel against the other dominant terms in the coupling c aa h to keep it small (the term proportional 
to A\ is multiplied by the singlet component of h or the doublet component of a, both of which are 
small in this region and hence A\ must be several TeV for this term to be significant). Exacerbated 
by the large A, these large values for A\ have a similar impact as large stop masses on the fine-tuning 
of the Higgs VEV and are therefore undesirable from the perspective of naturalness. However, the 
lower end of this range with |A^| ~ 1 TeV suggests that this scenario may be still compatible with 
the optimally natural region. 

The points with a bino-dominated LSP do not require large A\ and lie directly in the region of 
the NMSSM parameter space with minimal fine-tuning. Note that we find points with A spanning 
the full range over which we scan (see Table 3), which includes points both below and within 
the large-A regime (where the renormalization group running of A reaches a Landau pole before 
the GUT scale). In particular, all of these have small |/r e ff| ~ 200 GeV that both maximises the 
couplings Cax \2 an d c hxx 2 (through maximising the Higgsino fraction of X 2 within experimental 
constraints) and keeps the X 2 light. Note that in the bino-dominated case, c a ^ 2 and cannot 

be much larger than 0.2 because, unlike the singlino-dominated case above, the dominant terms in 
these couplings are proportional to g\ rather A. It is thus preferred from (3.20) that x.2 be as light 
as possible. Substantial contributions to the annihilation cross-section also arise from the s-channel 
annihilation via the second pseudoscalar 02 , which has mass typically in the range ~ 300 — 500 GeV 
because of the small values of p e g. The consequences of direct detection on the coupling c a xx for 
a singlet-like a described in Sec. 3.3 prevent the s-channel annihilation via the a from dominantly 
contributing to the cross-section (in contrast to the singlino-dominated case, where the /-c-dependent 
term in (3.14) allows c a ^ to be much larger), although we do not claim that our scan is exhaustive 
enough to discount this process altogether. Despite being more appealing, the bino-case is also more 
severely constrained, particularly by electroweak precision measurements and direct detection. 

In Fig. 4b, we show a scatter plot of the cross section for spin-independent DM-nucleon scat¬ 
tering versus the DM mass. We also plot the current limits from LUX and the projected limits 
from LUX with 300 live days of measurements [87], XENON1T [ 88 ] and LZ [89] (XENON10T 
has a similar projected sensitivity [ 88 ]) for a local DM density pg = 0.3GeV/cm 3 . Note that the 
uncertainties in the DM distribution invoked above may be correlated with the bounds from direct 
detection. Indeed, the range A £ [0.17,5.3] that we consider for the astrophysical uncertainty 
factor is determined by the uncertainties in the slope parameter 7 , the scale radius R s and the 
local DM density p & [11; 12], To illustrate how the direct-detection bounds are affected by the 
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71a [GeV] 

A k [GeV] 

M [GeV] 

A^eff 

l/^eff | 

Mi [GeV] 

A 

K 

tan /3 

range [-1000,1000] 

[-300,300] 

[300,600] 

±1 

[30,50] 

[0.8,1.4] 

[0.8,1.4] 

[1,4] 

step 10 

1 

100 


1 

0.05 

0.05 

1 


Table 4: Parameter ranges and step sizes used in the grid scan for the bino-dominated LSP anni¬ 
hilating into 66. For p, e g, the scan is over both possible signs. 


latter uncertainty, we plot them as shaded bands for local DM densities between p@ = 0.2 GeV/cm 3 
and p@ = 0.6 GeV/cm 3 . The LUX collaboration plans to update their analysis to 300 live days of 
measurements in the near future, extending the sensitivity by a factor of 5 [87]. This would test a 
considerable portion of the singlino-like points and almost the entirety of the bino-like points from 
our scan of the ha-channel. The vast majority of points with a singlino-like LSP would then be 
probed by XENON1T, although we find a few such points that would evade even the projected LZ 
limits. 

4.2 Parameter scan for the 66-channel 

We have performed a grid scan for the 66-channel and a bino-dominated LSP using NMSSMTools 
4.4.0 [32-34], We have fixed the masses of the squarks, sleptons, gluino and wino to the values 
mentioned previously and have scanned over the remaining free parameters of the scale-invariant 
NMSSM with ranges and step sizes as summarised in Table 4. The range for Mi is chosen such that 
the LSP is dominantly bino-like. We show scatter plots of the results from this scan in Figs. 5 to 7. 
As expected, we find that annihilation into 66 is the dominant channel. We demand that at least 80% 
of the annihilations are into 66 states to remove the very few points with a slightly smaller branching 
ratio to 66. The cross sections are restricted to lie in the range (a bb v)o £E [0.5,5] x 10~ 26 cm 3 /s 
required to explain the GCE discussed in Sec. 2. Green points have an LSP abundance in the range 
0.107 < IL^Ti 2 < 0.131 consistent with the observed DM abundance as discussed above. Red and 
blue points, on the other hand, have smaller relic densities. We keep the latter points because, in 
Sec. 5, we present a model in which NLSP decays can regenerate an under-produced LSP. We do 
not include points with abundances above the aforementioned range (which could become viable 
e.g. if there is significant entropy production after freeze-out). Blue points have a pseudoscalar 
mediator with a greater doublet fraction than singlet fraction, while for the red points the singlet 
fraction is larger. 

We have seen in Sec. 3.3.1 that constraints on the Higgsino fraction make it difficult to obtain 
the annihilation cross section required for the GCE without an additional enhancement. Such an 
enhancement arises if m a — 2so that the annihilation happens close to resonance. Focusing 
on this enhancement, we have scanned over the region 1/1.1 < 2 m^/m a < 1/0.9, where the ratio 
2 m^/m a measures the proximity to resonance. We show a scatter plot of this ratio versus the LSP 
mass in Fig. 5. Notice that only very few points come close to the upper or lower limit of our scan 
region (which approximately coincides with the plotted range) and that the vast majority of points 
lie closer to resonance, illustrating the necessity of this enhancement.' Points in Fig. 5a have cross 

7 It should be noted that higher-order quantum corrections to the pseudoscalar mass than those included in 
NMSSMTools 4.4.0 may perturb it sufficiently to significantly change the cross-section in the vicinity of the resonance. 
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Figure 5: Scatter plots of the LSP mass versus the ratio 2 m^/m a . Green points have an 
LSP relic abundance within 10% of the DM density determined by Planck, whereas for red and 
blue points it is smaller. Red and blue points have a pseudoscalar a which consists dominantly 
of sj and Hj , respectively. For the left panel (a), the late-time cross section is within the range 
{ a bb v )o £ [0.5,5] X 10~ 26 cm 3 /s, where the upper limit satisfies the published dwarf-spheroidal 
bounds from Fermi-LAT [50]. For the right panel (b), the upper limit satisfies the latest, preliminary 
dwarf-spheroidal bounds from Fermi-LAT [51], ( cr bb v)o £ [0.5,1.5] X 10 _26 cm 3 /s. 


sections in the range (a bb v)o £ [0.5,5] x 10 _26 cm 3 /s, where the upper limit is determined by the 
published dwarf-spheroidal bounds from Fermi-LAT [50]. Points in Fig. 5b are restricted to lie 
in the smaller range (cr bb v) o £ [0.5,1.5] X 10 _26 cm 3 /s, where the upper limit satisfies the latest, 
preliminary dwarf-spheroidal bounds from Fermi-LAT [51]. Note that we do not apply this second 
range to other plots in this section. 

A scatter plot of the cross section at freeze-out versus that at the present time is shown in Fig. 6. 
Notice that for many points, the former is orders of magnitude larger than the latter. This can be 
understood by noting that when m a is close to 2 m^, the annihilation is enhanced by the s-channel 
propagator at small velocities. Then if m a > 2m^, the annihilation can happen even closer to 
resonance during freeze-out due to the larger kinetic energy of the LSP in the thermal plasma. This 
enhances the cross section at freeze-out relative to that at the present time [90] and thus reduces 
the relic abundance (as follows from Fig. 6, these points typically have (cr bb v)T F ~ 0(1O 3 ) x (a bb v)o, 
leading to LSP relic densities fi^/i 2 < 10~ 3 ). Indeed, the majority of points in Fig. 5 that correspond 
to underabundant LSPs (in red and blue) have m a > 2 m^. Notice also that points with a singlet- 
dominated pseudoscalar (in red) tend to be closer to resonance than those for which it is dominantly 
doublet-like (in blue). This can be understood from the fact that the singlet does not couple to the 
SM and the doublet couples more strongly to neutralinos, so the product of couplings c a ^ • c abb 
becomes smaller with growing singlet-admixture to the pseudoscalar (cf. Eqs. (3.12) and (3.15)). 
On the other hand, notice from Fig. 5 that all points for which the LSP abundance accounts for 

However, we do not expect that these will qualitatively change our results, but would rather exchange points with 
otherwise similar characteristics in and out of the region satisfying the fits to the GCE. 
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Figure 6: Scatter plot of the freeze-out cross section ( a bb v)T F versus the late-time cross section 
(a bb v)o. Green points have an LSP relic abundance within 10% of the DM density determined by 
Planck, whereas for red and blue points it is smaller. Red and blue points have a pseudoscalar a 
which consists dominantly of s/ and Hj, respectively. 


the DM (in green) have m a < 2 m^. The larger kinetic energies at freeze-out then reduce the cross 
section for the process XX —> cl* —>• bb compared to the corresponding cross section in our galaxy. 
For the green points with LSP masses close to mz/ 2, this is counterbalanced by the resonant Z- 
boson mediated annihilation. This process does not contribute significantly to the cross section in 
our galaxy as it is helicity-suppressed in the s-wave and the width of the Z-boson is too broad, 
but becomes important in the early universe (where the p -wave becomes important) and raises the 
cross section to the required value to account for the DM. This effect was discussed previously in 
Refs. [17; 91; 92], 

Away from the Z resonance, both the cross section at freeze-out and today are simply dominated 
by a-mediation to 66-pairs, with the former being smaller than the latter. The freeze-out cross sec¬ 
tion is determined by the requirement that the LSP accounts for the DM, {<jv)t f ~ 2 x 10 26 cm 3 /s. 
Correspondingly, these points tend to have late-time cross sections that are larger than this, 
(a bb v)o > 3 x 10~ 26 cm 3 /s. These large cross sections are still compatible with the GCE because 
of the relatively large uncertainties in the fits in Tables 1 and 2 and the large astrophysical un¬ 
certainty factor A. Indeed, if we impose the preliminary dwarf-spheroidal bounds as in Fig. 5b, 
(a bb v)o < 1.5 x 10~ 26 cm 3 /s for m DM ~ 30 — 50GeV [51], only green points with masses close 
to mz /2 survive. Similarly, points for which the LSP is underabundant (in red and blue) with 
m a < 2 niy. in Fig. 5 have relatively large late-time cross sections. For these points, however, the 
thermal broadening of the resonance does not sufficiently reduce the freeze-out cross section to 
account for the DM. 

In Fig. 7a, we show a scatter plot of the cross-section for spin-independent LSP-nucleon scatter¬ 
ing versus the LSP mass. As discussed in Sec. 4.1 for Fig. 4b, we plot the current limits from LUX 
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Figure 7: The left panel (a) shows the spin-independent cross section for scattering of LSPs on 
nucelons versus the LSP mass, while the right panel (b) shows the branching ratio of invisible Higgs 
decays versus the LSP mass. Green points have an LSP relic abundance within 10% of the DM 
density determined by Planck, whereas for red and blue points it is smaller. Red and blue points 
have a pseudoscalar a which consists dominantly of sj and Hj , respectively. 


and the projected limits from LUX with 300 live days of measurements [87] and XENONIT [88] as 
shaded bands covering local DM densities between p Q = 0.2 GeV/cm 3 and p 0 = 0.6 GeV/cm 3 . The 
LUX collaboration will test a considerable portion of the points from our scan for the 66-channel 
in the near future. The remaining points will be tested by XENONIT. 

In Fig. 7b, we show a scatter plot of the branching ratio for Higgs-to-invisible decays versus 
the LSP mass. The better precision in the measurement of the Higgs properties projected for 
Run 2 of the LHC and after the high-luminosity upgrade will allow for another way of probing 
our points. Indeed, the CMS collaboration expects to probe branching ratios for invisible Higgs 
decays at the level of Br (h —> inv) E [8, 11]% at lcr for an integrated luminosity of 300 flu 1 and 
Br(/t —> inv) E [4,7]% at lcr for 3000 fb _1 [93]. We have plotted these ranges as colored bands in 
Fig. 7b. As one can see, if no significant deviations from a SM-like Higgs are found at the LHC 
during the high-luminosity run, only a small slice of parameter space with DM masses in the range 
46 GeV < < m^/2 would survive. 

To conclude, we find a small number of points (around 5 % of all points which do not overdose 
the universe) that can simultaneously accommodate the GCE via the 66-channel and the observed 
DM abundance from thermal freeze-out. The majority of our points, however, have an underabun- 
dant LSP. In the next section, we discuss a model where the right abundance of LSPs is obtained 
from late NLSP decays. This introduces a new degree of freedom that would allow for the observed 
DM abundance to be separately accounted for from the GCE. The majority of the viable points 
will be tested by LUX [89] this year, while the remainder will be probed by XENONIT [88]. 
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5 Seesaw extension with sneutrino LSPs or NLSPs 


5.1 Seesaw in the scale-invariant NMSSM 

We consider a simple extension of the NMSSM, the seesaw [94-98], that provides new possibilities 
for the dark sector and simultaneously addresses the generation of neutrino masses absent in the 
original model. We introduce three right-handed neutrinos and their scalar superpartners, i.e. 
gauge-singlet superfields W (the i denotes the flavour index). The NMSSM superpotential is 
extended to include the terms 


(5.1) 


Wseesaw — A jv SNN + y N L ■ H u N , 


seesaw 


and the soft supersymmetry breaking terms are 



where we have suppressed the flavour indices on the supermultiplets N and L, along with those in 


the parameters Ajv, y N , m \i A\ N and A Vn . 

The Majorana masses mj\r = 2Xnv s of the fermionic components are generated by the YEV of 


the scalar component of the singlet S and because of this link with the Higgs sector are naturally 
at most of order TeV. Provided that the Majorana masses are hierarchically larger than the Dirac 
masses, the seesaw mechanism yields active neutrino masses 



(5.3) 


In order for the active neutrino masses to be smaller than ~ eV, the Yukawa couplings then need 
to be \y N \ < 1(T 6 . 


As y N is much smaller than the other couplings, the model is approximately invariant under 
a Z 2 parity under which only the right-handed neutrino supermultiplet is odd (cf. Eqs. (5.1) and 
(5.2)). Processes that involve odd numbers of these particles are therefore suppressed by the small 


coupling y N , whether mediated by the Yukawa coupling or mass mixing. If the lightest of the right- 
handed sneutrinos is not the LSP, this approximate symmetry can have the interesting implication 
that this particle is relatively long-lived. We will later exploit this property in connection with the 
non-thermal generation of DM. 


In order to simplify the discussion, we restrict to just one flavour of right-handed neutrinos in 


the following. We assume that CP is conserved in the sneutrino sector and take all couplings and 


masses to be real. The real and imaginary parts of the sneutrino N = (Nr + iNi)/V 2 then do 
not mix with each other. Mass mixing with the left-handed sneutrino sector is proportional to the 
small coupling y N and can be ignored. The masses of the right-handed sneutrinos are then given 
by 


= m 2 ^ + 2X n v 2 ± Xn (An 2 sin2/3 - 2A\ N v s - 2nv 2 ), (5.4) 


(5.4) 


where the + and — correspond to Nj and Nr respectively. These masses can be regarded as 


independent free parameters, simply exchanged for the parameters m and A\ N in the Lagrangian. 
This in particular fixes 
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(5.5) 





Figure 8 : Feynman diagram contributing to 
bosons. 



sneutrino annihilation to SM fermions and gauge 


The size of A\ N is then partly controlled by the mass-splitting between the CP-even and CP-odd 
sneutrino. 

The possibility that the sneutrino is the LSP has been discussed previously in [35] and [36] , the 
latter of which also provides a cursory illustration of compatibility with the GCE. Because scalar 
DM respects Bose-Einstein statistics and its annihilation is not p -wave suppressed, it is easier to find 
amenable regions of parameter space than for the neutralino, although direct-detection constraints 
are still strong. Fits to the spectrum of the GCE for annihilating sneutrino DM (including scalar 
annihilation channels) were examined more closely in [37]. Collider signatures of this scenario were 
discussed in [99]. 

In Sec. 5.2, we begin by exhibiting the case of sneutrino DM annihilating dominantly into Higgs 
scalars and pseudoscalars. Then in Sec. 5.3, under the assumption of an especially small coupling 
y N , we make use of the approximate Z 2 invariance in the right-handed sneutrino sector to point-out 
the distinct possibility of metastable sneutrino NLSPs that decay into neutralino LSPs. This is 
applicable to the parameter regions in Sec. 3 where the annihilation cross-section of the neutralino 
LSP in the early universe is too large. 

5.2 Sneutrino LSP 

We now analyse the possibility of sneutrino DM to explain the GCE. We focus on annihilation 
into Higgs scalars and pseudoscalars, with the mass and cross-sections prescribed in Sec. 2 , com¬ 
plementing the recent analysis for smaller values of A and k in Ref. [37]. The case of annihilation 
into light quarks has been discussed previously in Ref. [36]. For simplicity, we consider only one 
flavour of right-handed (s)neutrinos and neglect flavour mixing. Furthermore, we assume that the 
mass eigenstate Nj is the LSP (henceforth denoted N except where it is necessary to distinguish it 
from the other sneutrino), although the subsequent analysis would be similar if Nr was the LSP. 

The couplings of the right-handed (s)neutrino to the Higgs, the lightest pseudoscalar and the 
neutralinos are listed in Appendix A. Note that these couplings are controlled by either A at or 
A nA\ n . The coupling Atv is constrained in several ways, the first of which we discuss is by direct- 
detection experiments. Analogously to Sec. 3.2, we assume isospin symmetry and obtain [75; 100] 

m u <5l2 c hNN m d ‘-’ll C hNN Ir. 

a u = - 7 =-- — i , ad = - 7 =- —— i (5-6) 

2 V 2v sin p 2m^m^ 2y/2v cos p 2mf l m^ 

for the coefficients which enter in Eqs. (3.5) and (3.6). Assuming a Higgs with SM-like couplings and 
that no other scalars contribute to LSP-nucleon scattering, for a sneutrino mass ~ 100 GeV, 
the bound on the direct-detection cross section by LUX gives \n < 0.1. 
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The annihilation of sneutrinos into SM fermions and gauge bosons is mediated via s-channel 
Higgs exchange, as shown in Fig. 8 . The cross section is given by 


(°sm v )o 


_ C hNN _ 

m jv( 4m ^ - m h ) 2 


Tfe( 2 mjv), 


(5.7) 


where c h f fl y is the coupling of two sneutrinos to the Higgs and Yh{2m^) is the decay width of 
the SM Higgs evaluated for a Higgs mass of 2 m^. For > rri\y , the dominant process here is 
NN —> h* —> W + W~ , due to the large gauge couplings. The cross section is controlled by the 
coupling which is given in terms of the couplings and soft masses in Eq. (A.6). It scales 

with the dimensionful parameters v, v s and A\ N . The latter two are multiplied by the singlet 
fraction of the Higgs which needs to be small in order to be compatible with the measured Higgs 
couplings, typically 5i3 ~ 0.1. Since we restrict ourselves to v s and A\ N below TeV, the contribution 
proportional to v dominates the coupling and we can approximate it as 

\ c hNN\ ~ A at A v sin 2/3 « 1-5 GeV ^ ^ (5.8) 

for a SM-like Higgs with Sn ~ cos /3 and S 12 ~ sin /3. Demanding that the cross section for 
annihilation into SM particles be much smaller than the one required to explain the GCE, ((j gM n)o <C 
2 • 10 -26 cm 3 /s, gives the bound | c /j,tvtv I < 20 GeV for an LSP mass above the Higgs resonance and 
up to ~ 10 GeV below the W threshold. The bound is even stronger above the W threshold, with 
I c ft .jvtv I 1$ 5 — 13 GeV for masses = 100 — 150 GeV. This translates into the bound A n < 0.1. 
Note that this bound becomes stronger if the astrophysical-uncertainty factor A is large and smaller 
cross sections are sufficient to generate the GCE (see Sec. 2). Also note that, for < m\y, the 
largest branching fraction is into 66-pairs, although the branching fraction into VFIW-pairs is not 
much less important except for a small range of masses above the Higgs resonance that we do not 
consider further (below which constraints from invisible Higgs decays into sneutrinos give a bound 
on c h ^fj that is too strong for (5.7) to give a large enough cross section for the GCE). 

The dominant processes that allow for the annihilation of sneutrinos into scalars or pseudoscalars 
from the Higgs and singlet sector are 1) via the four-point scalar coupling, 2) via s-channel mediation 
of the Higgs and 3) via t- and u-exchange of a sneutrino. We show the Feynman diagrams for these 
processes in Fig. 9. The cross section for the annihilation of sneutrinos into Higgs pairs is given by 


(Phhv) 0 


\l m2 N- m2 h c ftjVJV c hhh 

647r m 3 ~ hhNN 4 m 7 — rnj 

N N h 


9 c 2 2 

wC hNN 

ml — 2m 2 - 

h AT 


(5.9) 


where c hh ^jy i s the coupling between two Higgs and two sneutrinos and Chhh is the triple Higgs 
coupling. We find that it is difficult to obtain the required cross section to generate the GCE via 
the hh -channel while keeping annihilation into SM fermions and gauge bosons subdominant. This 
is due to the fact that the s-, t- and n-channel processes are determined by the same coupling 
as the annihilation into SM particles, while the coupling c hh ^ 1 y relevant for the direct annihilation 
is also controlled by Aat < 0.1 (see Eq. (A. 8 )). In addition, the annihilation into a Higgs pair is 
phase-space suppressed because it needs to happen close to threshold for a good fit to the GCE. This 
makes it difficult to enhance the 6 ,/ 1 -channel over annihilation to SM particles. We have performed 
a short uniformly random scan of the parameter space in the natural region of the NMSSM with the 
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Figure 9: Feynman diagrams contributing to sneutrino annihilation to hh- or aa-pairs. Note that 
there is a u-channel diagram in addition to the f-charmel one. 


LSP mass constrained to lie within 30 GeV of the Higgs mass threshold. However, the branching 
fraction into hh was never found to exceed ~ 40% (with the remaining annihilations dominantly 
into W + W~). Although we did not search more extensively and cannot conclusively rule this 
possibility out, we will not consider it further. 

Annihilations into singlet-like scalars or pseudoscalars are more promising for generating the 
GCE because the sneutrino can couple more strongly to relatively pure singlets than to doublet-like 
states. For a SM-like Higgs with small singlet fraction, the largest terms in the couplings c h ^jy 
and c hh fjpj (see Eqs. (A.6) and (A. 8)) are proportional to A sin 2/3, whereas the terms proportional 
to the singlet fraction can make the corresponding couplings for a singlet-like state larger without 
affecting the annihilation cross-section into other SM particles such as W + W~. Scenarios with a 
light, CP-even singlet-like state (possibly lighter than the Higgs) were studied in [37] for small A 
and k < 0.1. However, such a light singlet is difficult to arrange in the more natural region of the 
NMSSM with large A and k as this state typically mixes with the Higgs and the mixing is enhanced 
with A, leading to large deviations in the Higgs couplings from SM values. The mixing must then 
be tuned small in order for the Higgs couplings to be consistent with experiment (see Eq. (4.1)). 

We will instead focus on the annihilation channel of light pseudoscalars. These can arise natu¬ 
rally as the pseudo-Nambu-Goldstone boson of an approximate, spontaneously broken P-symmetry 
in the Higgs sector if A\ and A K are relatively small. The thermally-averaged annihilation cross 
section for this process is 


{&aaV)T F 



64 7r m 3 ~ 


C aaNN 


C hNN C haa 

4 "4 - m h 


2 c 


2 _ 

aNjN R 


9 2 2 

mi — rn~ — m- 
a Nj n r 


2 


(5.10) 


where Ch aa is the Higgs coupling to two pseudoscalars and c a jr / n h , c aaNN an d c hNN are cou ~ 
plings of two sneutrinos to the pseudoscalar, two pseudoscalars and the Higgs, respectively. We 
have assumed that the cross-section is well-approximated by its velocity-independent component, 
although care must be taken in performing the thermal-average close to the annihilation threshold 
[90; 101] (we also ignore the Higgs resonance). 

We will now determine regions for the couplings \n and A\ N (or, equivalently, m^ R ) where 
(cr aaHjTp ~ 3 x 10 _26 cm 3 /s can be achieved so that the observed DM abundance is obtained from 
thermal freeze-out. We will use these regions for our scan below. To this end, we determine 
parameter values for which each term in (5.10) individually gives the required freeze-out cross- 
section. Since generically all three terms contribute, this gives approximate upper bounds for the 
corresponding parameters where this annihilation channel can be viable. The first two terms provide 
suitable regions for Atv, while the third points to an appropriate range for m^ R . Let us first consider 
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the case that the four-point scalar coupling dominates the cross section. The expression for c aa ^ 
in terms of the couplings and soft masses is given in (A. 10), where V is a diagonalisation matrix 
that relates the mass eigenstates at of the pseudoscalars to the interaction eigenstates according to 



V 



(5.11) 


Assuming a very singlet-dominated pseudoscalar with V\z ~ 1 and Vii,V \2 1, we find c aa ^^ ss 
— 2Xj\f(2X]\r + k). The cross section does not exceed that required to obtain the observed DM 
adundance when A at ;$ 0.05 for k « 1, « 100 GeV and m a ~ 90GeV. This bound is relatively 

independent of the values of and m a in the region of interest or becomes stronger (provided 

that the masses are not too close to the threshold, in which case (5.10) becomes invalid). This 
satisfies the bound above required for SM channels to be subdominant. A similar argument for the 
case where the s-channel process dominates the cross section gives a comparable limit on A at. 

Via c a fj i . the third term in (5.10) depends upon A nA\ n multiplied by the large singlet 
component of the a (see Eq. (A.ll)). Since we fix A nA\ n in terms of the sneutrino mass (see 
Eq. (5.5)), we can obtain a constraint on this mass from the requirement that the third term in 
(5.10) does not exceed the required freeze-out cross section. To this end, we take the limit Aat —> 0 
keeping Aa j A\ n fixed (in which case the interference from the other amplitudes is negligible) and 
assume that m 2 - 3> m 2 ~ ,m% so that the terms proportional to m,-, dominate. In this limit, 

Nr Ni a Jr" Jr- A/ fl 


C aNj Nr 


m 


Vn 


Ne 


y/Si 


and the sneutrino annihilation cross-section is 


( &aa,v)T F 



^13 m N R \ 4 


(5.12) 


(5.13) 


The cross section does not exceed that required to obtain the observed DM abundance when 
m N R < 0.75u s /|"Pi3| for m^ li ~ 150 GeV and m a ~ 0.9 m^. This is roughly independent of 
and m a within the best-fit range. We have v s < 500 GeV in the natural region of parameter space, 
so if \V n \ > 0.9, then < 400 GeV. We caution the reader again that this analysis is only 

approximate, but reasonably describes our target parameter space. 

We have performed a random scan of the parameter space of the scale-invariant NMSSM ex¬ 
tended by right-handed neutrinos. NMSSMTools 4.3.0 and micrOMEGAs 3.0 were used to calculate 
the spectrum, phenomenological constraints, DM signals and relic abundance. The code was mod¬ 
ified to calculate the relic abundance of the sneutrino LSP, using modified model files for CalcHEP 
[102] generated with LanHEP [103]. All calculations of neutrino/sneutrino cross sections were per¬ 
formed at tree-level. The interactions mediated by the coupling y N were omitted as they contribute 
negligibly for all processes of interest here. The ranges for the parameters used in the scan are 
given in Table 5. To concentrate on a light pseudoscalar, we have only scanned over relatively small 
values for A\ and A K . According to the heuristic analysis above, we have restricted \Xn\ < 0.1 
and m^ R < 400 GeV. To ensure that the LSP mass lies in the best-fit region of the aa-channel 
as in Fig. Id, m^ was restricted to uniformly random values between [m a ,m a + 30GeV], up to a 
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Ax [GeV] 

A k [GeV] 

Ws\ [GeV] 

A 

k tan/3 Am^, ( [GeV] 

Xn 

m N R [ GeV ] 

[-100,100] 

[-10,10] 

[200,500] 

[0.5,1.5] 

[0.2,1.4] [1.2,4] [0,30] 

[-0.1,0.!] 

o 

O 


Table 5: Parameter ranges for the random scans for the aa-channel, where A — m a . 
Both signs of fj, e g are scanned over. 



(a) (b) 

Figure 10: Results from the random scan with parameters given in Table 5. The left panel (a) 
shows the annihilation cross section at the present time, (cr aa v)o, and at freeze-out, {<jv)t f - The 
right panel (b) shows the spin-independent cross section for scattering of LSPs off nucleons versus 
the LSP mass. 


maximum of 150 GeV. The remaining NMSSM parameters were fixed as in Sec. 4 with the addition 
that M\ = 2 TeV. Likewise, collider and flavour constraints were also imposed as in Sec. 4. We 
also demand, as before, that the LSP has the necessary relic abundance to entirely account for 
DM, satisfies the constraints imposed by the LUX experiment and has an annihilation cross section 
consistent with the GCE. As we find A < 1.1 for all points presented and no clear preference for 
tan /3 > 3, we expect all points to be consistent with electroweak precision tests and have refrained 
from explicitly checking this (see also Ref. [67]). We emphasise that our approach is not an exhaus¬ 
tive scan of the parameter space, but rather an illustration of how this model provides a natural, 
supersymmetric account of the GCE. 

The results from our scan are presented in Fig. 10. All points have a branching fraction for 
annihilations into aa of at least 80%. As suggested from the above analysis, most viable points he 
in the range |Ajv| 1$ 0.04. The appeal of this scenario is that the cross-sections are not severely 
constrained by experiments and are not dependent on special effects, like a resonance. The latter 
point becomes clear in Fig. 10a, where we plot the cross-sections at freeze-out and late-times. These 
are of comparable size for the vast majority of points. We find that the few outlying points with 
large (a aa v)o <1 3 x 10~ 26 cm 3 /s have smaller cross-sections in the early universe because the total 
freeze-out cross section has additional contributions from coannihilation with Nr, while the points 
with larger (<jv)t f are resonantly enhanced by a singlet-like Higgs scalar. Just as for the case of the 
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neutralino annihilating into ha, the typical cross sections for most of the points are close to half of 
that prescribed for the best-fit regions in Fig. Id and therefore require an astrophysical-uncertainty 
factor A ~ 2 to explain the GCE. We have plotted our estimate for the limit on (a aa v )o from the 
latest Fermi-LAT study of dwarf spheroidals, that we have discussed in Sec. 2, as a dashed line in 
Fig. 10a. The vast majority of points satisfies this limit. We show a scatter plot of the cross section 
for spin-independent scattering of LSPs off nucleons versus the LSP mass in Fig. 10b. As discussed 
in Sec. 4.1 for Fig. 4b, we plot the current limits from LUX and the projected limits from LUX 
with 300 live days of measurements [87], XENON1T [88] and LZ [89] as shaded bands for local 
DM densities between = 0.2GeV/cm 3 and Pq = 0.6GeV/cm 3 . Notice that the vast majority of 
points would satisfy the projected, updated LUX bounds, whereas only very few would evade the 
LZ bounds. 

5.3 Metastable sneutrino NLSP 

As described in Sec. 4.2, the pseudoscalar resonance can enhance or suppress the annihilation cross- 
section of the neutralino LSP in the early universe compared to that at late times. For regions 
of parameter space that explain the GCE via the 66-channel and satisfy m a > 2 m^, the thermal 
relic abundance is then too small to account for the DM content of our universe. We now discuss 
the possibility that a metastable right-handed sneutrino NLSP freezes-out and decays into the 
neutralino LSP to produce the observed DM abundance. This realises a similar simplified DM 
scenario proposed in [104], albeit in a different context. 

We will again restrict our discussion to one flavour of right-handed neutrinos and will neglect 
flavour mixing. It will turn out that the suppression of flavour mixing is a necessary consequence of 
the requirement of metastability of the sneutrino. Furthermore, we will continue to consider only 
the CP-odd mass eigenstate of the sneutrino. Another qualitatively similar possibility that we will 
not explore is that the CP-even or both states are metastable. 

The right-handed sneutrino NLSP freezes-out before the LSP does because of the smaller LSP 
mass. In order for the NLSP to contribute to the LSP relic abundance, it needs to decay after 
the LSP freezes-out. In addition, the NLSP lifetime must be sufficiently long to ensure that LSP 
annihilations do not resume once these are repopulated (the “freeze-out and decay” scenario in 
[105]). Neglecting changes in entropy due to the late sneutrino decays, the total density fraction of 
DM is then given by 

TTl ~ 

Ddm = -I-(5-14) 

Here Ddm is the DM density fraction, is the density fraction of LSPs produced as a thermal 
relic and is that of NLSPs expected today had they not decayed. For the points of interest 
from Sec. 4.2, ^<0^. 

The requirement that the sneutrino be long-lived places an upper bound on the coupling y N that 
determines the non-conservation of the approximate Z 2 in the right-handed neutrino sector. We 
will investigate this bound below. However, Z 2 -preserving decays into the right-handed neutrino, 
such as IV —> x + -/V, are controlled by the coupling Ajv- This cannot be much smaller than 
the electroweak gauge couplings for a sneutrino of mass ~ GeV — TeV to freeze-out with a relic 
abundance that reproduces the observed amount of DM. 8 These fast decays must be kinematically 

8 We do not consider alternative non-thermal production possibilities that involve a very small Aat such that the 
sneutrino is never in thermal equilibrium with the plasma to begin with, see [106] for an example in the MSSM. 
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forbidden if the sneutrino is to be long-lived and therefore the sneutrino mass must satisfy 


rrifi < + m n 


m x + 


2 AAr/igff 

A 


(5.15) 


The neutralino mass is ~ 30 — 50 GeV in the parameter region of interest for the 66-channel 
(see Sec. 4.2). As we focus on A ~ 1 to improve naturalness, we therefore require that < g e s, 
taking Aat ~ 0(1), where g e g is typically of order a few hundred GeV. The dominant sneutrino 
decay mode then is N —> x + v, which has a rate 


T(iV-» x + v) = 


~Nxv\ 

8n 


m N 1 - 


m 


m 


x. 

~N, 


where 


c Nxv ~ ~ 2 9lJ ^ uNP 


The lightest sneutrino N has a left-handed sneutrino component given by 

Vi v 


N v « (2\nv v s sin f3 + v(Xv s cos /3 — A y sin /3)) 


2 2 ’ 
- m u 


(5.16) 

(5.17) 

(5.18) 


where ma is the mass of the left-handed sneutrino. The mixing factor given by the second term in 
(5.17) is therefore also directly proportional to y N <C 1. 

To ensure that the replenished LSPs remain frozen out of chemical equilibrium, the Hubble 
rate must be faster than the annihilation rate of the LSPs as the sneutrinos decay. Calling the 
temperature at which the sneutrino decays with its mean lifetime, this requires that 


( av ) Td Yoo s(T d ) < H(T d ). (5.19) 

Here, (av)T d is the thermally-averaged annihilation cross section of the LSPs evaluated at tempera¬ 
ture Td, Yoo = n*(T 0 )/so where n^(To) is the present-day LSP number density and so = 2.22 x 10" 38 
GeV 3 the present-day entropy density of the universe and H(Td ) = irx/g*/90T%/Mpi and s(Td ) = 
(27r 2 /45)g*Tj are respectively the Hubble rate and entropy density at the time of the decays with 
g* the active relativistic degrees of freedom. The time interval between sneutrino freeze-out and 
decay could encompass the QCD phase transition, over which the number of degrees of freedom in 
the plasma must be carefully treated. We include this temperature dependence in our analysis, but 
for simplicity will treat this as a constant g * in the analysis presented below, which is sufficient for 
an order of magnitude estimate. 9 

The inequality (5.19) provides an upper bound on the decay temperature of the sneutrinos: 

T d < \/ g 4 9 5 v x • (5-20) 

V 87 x z g* Yoo A/pi {crv)T d 

9 Note that the case of equality in (5.19) approximately describes the special possibility of “reannihilation” [105]. 
In this scenario, as the sneutrinos decay, the neutralino number density grows sufficiently for annihilations to resume. 
These annihilations continue rapidly until the decays finish, after which the neutralinos quickly freeze-out again when 
their number density falls to the critical level given by the equality in (5.19) (technically the temperature of LSP 
freeze-out is actually slightly lower than that corresponding to the sneutrino lifetime). The resulting neutralino relic 
abundance is approximately independent of the original sneutrino number density, provided that the latter is at least 
as large, because the surplus neutralinos rapidly reannihilate. 


33 









Note, however, that (av)T d is itself a function of the temperature Tj. Indeed, for the points in 
parameter space that we consider here, the pseudoscalar resonance results in an enhancement of 
the annihilation cross section during freeze-out compared to late-times. The cross section decreases 
at temperatures below the LSP freeze-out temperature. Using parallel analysis to that of [101], the 
cross section may be well-approximated in the vicinity of the resonance as 

Mt ^ 1 (^) 1 erfc( - i^m^z R /T)) . (5.21) 

Here zr = 6r + i^R with er = (mj — 4m|)/4m~ and 7 r = m a r a /4m~, and T a is the decay width 
of the pseudoscalar. Using this, comparison of (5.19) at different temperatures reveals that it is the 
closing of the thermal-broadening of the resonance with the cooling of the plasma that determines 
the temperature after which the neutralinos will remain frozen-out. From Sec. 4.2, typical values 
are er ~ 0.015, T a ~ IMeV, |c ab 5 | ~ 0.05 and |c a ^| ~ 0.01 (where the cross section is most 
sensitive to er and r o ). Equation (5.20) then reduces to Tj, < 0.1 GeV. This bound can weaken by 
up to an order of magnitude if the masses are further from resonance, but is otherwise relatively 
insensitive to variations in the parameters over their ranges considered in our scan. 

This bound on the decay temperature Td corresponds to a bound on the decay rate T^ of the 
sneutrino. Denoting the time and temperature of sneutrino freeze-out respectively by F and 
T* F and the time of sneutrino decay by td, we have 


1 _ + + _ fQOMpi ( 1 1 ^ 

> ~ “ * F “ V* - ( t ? Ti J- 


(5.22) 


With Tjj F ~ rrijj /25 ~ few GeV, the bound Td < 0.1 GeV gives < 10 20 GeV. Using (5.16) 
and neglecting the second term, this translates into y N < 1CU 8 for rrij 7 ~ 1.1 • m ^ and |iVi 4 | ~ 0.01 
(as we typically find for a bino-like LSP, see Sec. 3.2), with a stronger bound for larger sneutrino 
mass. To ensure that the decays occur before big bang nucleosynthesis, we require that > 1 MeV. 
This corresponds to T^ > 10 ~ 25 GeV and y N > 1CT 12 . 10 

Recall that we have restricted our discussion to one flavour of right-handed sneutrinos. In 
general, the sneutrino may decay into any of the active neutrino flavours, in which case the decay 
width (5.16) and coupling (5.17) generalise and the bound on T^ derived above applies to the 
sum of partial decay widths for each flavour. The bound on the coupling c^ v which follows from 
(5.22) then applies to the sum of each of these generalised couplings added in quadrature. The 
resulting bound on the coupling y N derived above then applies to each of the Yukawa couplings of 
the superfield containing the NLSP to the active neutrino flavours. This implies that the neutrino 
superpartner of the NLSP cannot significantly contribute to the mass generation of the active 
neutrinos. The other two flavours then need to have corresponding couplings y N ~ 1CU 6 in order to 
reproduce the known mass splittings among the active neutrinos. Furthermore, this implies that 
flavour mixing between the NLSP and the other two right-handed sneutrino flavours must also be 
sufficiently suppressed to prevent premature decays. Note that this includes the suppression of 
flavour mixing that originates from the couplings Aat and A\ N and the masses rri'A. This means 
that the right-handed neutrino superfield that contains the NLSP is effectively decoupled from the 

10 As this is peripheral to our analysis, we do not carefully consider the possibility that big bang nucleosynthesis 
may be compatible with an epoch of decaying sneutrinos. 
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other two flavours and the lepton doublets, thereby justifying in hindsight our restriction to one 
flavour. The Z 2 parity mentioned in Sec. 5.1 may be restricted to this decoupled flavour. 

The freeze-out and decay of the N may be complicated by its interactions with the N and the 
X- Furthermore, the N may also be metastable and freeze-out with the N and x, making this an 
effectively multi-component DM model. Note that the Nr can similarly complicate the analysis if 
it is metastable. However, here we will assume that its mass is sufficiently high so that it decays 
on short timescales and can be neglected. The behaviour of the N must be carefully treated. If 
the mass of the N satisfies mjv > then it decays rapidly into the N and x and plays 

no role in DM production. This occurs for large fx e g or Aa j (for example, for the parameter values 
used in the scan in Fig. 11, the right-handed neutrino is short-lived and the sneutrino freezes-out 
with the necessary abundance for A jy « —0.15 and ~ 69 GeV). However, if these fast decays 
are kinematically forbidden, the right-handed neutrino will also be metastable and have a lifetime 
similarly determined by the small coupling y N . This can complicate the production of DM in 
several ways. 

Firstly, the self-annihilation cross-section of the right-handed neutrino is typically smaller than 
that of the sneutrino by one or two orders of magnitude. The neutrino will therefore freeze-out 
with a number density greater than the sneutrinos by roughly the same factor (assuming that the 
s-wave approximation holds). If kinematically permitted, a non-negligible fraction of these can 
decay into pairs of neutralinos through, for example, the decay chain N —>• va followed by a —> XX , 
thereby contributing considerably to the LSP density. A possibility that we leave for future work 
is that the LSPs are entirely produced by the decaying neutrinos while the sneutrinos freeze-out 
with negligible density. 

Secondly, the N and the N could affect their respective relic densities through “semi-annihilation” 
processes such as NN —> x a that could keep them in equilibrium for longer. 11 If these processes 
are significant, the departure of the N from chemical equilibrium during the freeze-out of the N 
would inextricably couple the dynamics of the thermal production of both species. In this case, a 
reliable calculation of the relic abundance of the sneutrino would require solving a pair of coupled 
Boltzmann equations. This is beyond the scope of this paper. Let us, however, point-out that the 
relic abundance of the N can be reliably calculated from the self-annihilation rate alone if the N 
freezes-out at a lower temperature and the “semi-annihilation” rate of the N with the N (and, for 
that matter, with the x) is subdominant. Although the N and x are lighter and hence much more 
abundant than the N, the semi-annihilation rates may be nevertheless mildly suppressed compared 
to the self-annihilation rate of the sneutrinos. This is both because the cross sections are smaller 
(typically suppressed relative to the self-annihilation cross-section by the small singlino fraction of 
the neutralino) and the metastability condition (5.15) in conjunction with the condition m a ~ 2 
implies that the reaction is often close to or below threshold. 

We have performed a scan in the sneutrino parameter space for an example point in the NMSSM 
for which a bino-like LSP could be responsible for the GCE but would otherwise be thermally 
underproduced (see Sec. 4.2). We have used NMSSMTools 4.3.0 and micrOMEGAS 3.0 with the 
same set-up as in Sec. 5.2 to calculate the relic abundance of sneutrinos under the assumption 
that these particles are stable (again neglecting y v -mediated interactions and assuming that these 

11 The destruction of a single N must either consume or produce a right-handed neutrino or be suppressed by 
the small coupling y N . Due to It-parity conservation, the y or a heavier 7?-odd particle must also be produced or 
consumed in such a reaction. Hence our use of terminology - the N, N and x ma y be regarded as particles in a 
Z 2 x Z 2 dark sector. 
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Figure 11: Scan over A n and m^ for an LSP of mass 63 GeV. The NMSSM parameters are selected 
as tan (3 = 3, A = 0.8, k = 1, fi e g = —350 GeV, A\ = —10 GeV, A K = 12.86 GeV and M\ = 60 GeV, 
which give m a = 125.3 GeV. The remaining sparticle and gaugino masses were fixed as in Sec. 4.2 
and we set = 300 GeV. This region in parameter space has a bino-like LSP which can 

explain the GCE but would be underabundant in the absence of the sneutrino NLSP. Dark grey 
denotes regions where the sneutrino freezes-out with a relic density that produces x DM consistent 
with the density determined by Planck. The sneutrinos are underproduced in the white regions 
and overproduced in the inner light grey regions. The shaded wedge violates the metastability 
condition, while the invisible branching fraction of the Higgs exceeds the experimental bound in 
the shaded slice in the upper-left corner. 


are small enough to satisfy the criterion above). However, this calculation incorporates only self- 
annihilations in the Boltzmann equations (and coannihilation with heavier i?-odd particles) and 
does not include the possibility of semi-annihilations. To check for the latter, we have compared the 
semi-annihilation rates for the reactions NN —> \a (or h) and N\ — > Na (or h) (assuming that the 
N is still in chemical equilibrium during the freeze-out of the N ) to the self-annihilation rate of the 
N at the freeze-out temperature computed by micrOMEGAS 3.0, using cross sections calculated by 
CalcHEP that were then thermally averaged. Likewise, the N self-annihilation rate was computed 
and compared to verify the assumption that it was in equilibrium at the beginning of freeze-out 
of the N. Where this holds and the semi-annihilation rates are smaller than the sneutrino self- 
annihilation rate by at least a factor of 10, we can safely neglect the impact of the N and x on 
the freeze-out of the N (the number densities of the lighter N and x will decrease faster than that 
of the N whether in equilibrium or freezing-out, so the semi-annihilation rates will also diminish 
faster as the temperature falls). We provide an example region where this occurs in Fig. 11 as a 
proof of principle. Over this region, the N is metastable and freezes-out after the N, while the 
semi-annihilations are weak enough to be ignored. Furthermore, its mass satisfies mjv < 2so 
that decays into the LSP are kinematically forbidden. The mass of the Nr, on the other hand, is 
large enough to ensure that it is short-lived by Ajv-mediated decays. Invisible Higgs decays into 
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sneutrinos or right-handed neutrinos can provide a further constraint on the parameter space if 
these species are light enough. As before, we require that the branching fraction for invisible Higgs 
decays satisfies Br(/r —> inv) < 0.24 [76; 77], which we compute at tree-level. 

We also check the possibility of conversion of the N into a x by Z 2 -violating annihilation with 
active neutrinos v mediated by the small coupling y N . Despite this small coupling, the much greater 
abundance of the v could enhance the rate for this process relative to the self-annihilation rate of 
the N. The largest of the processes mediated by the couplings y N is Nv —> %a, where the a 
may be on- or off-shell depending upon the total centre-of-mass energy. However, a rough estimate 
indicates that this process can be safely neglected. Indeed, the sneutrino number density in the non- 
relativistic regime in equilibrium is n^ ~ (rn^T )2 exp(—^A) and the left-handed neutrino number 
density is n v ~ T 3 . Ignoring this process, freeze-out occurs at temperature m^/T ~ 20 — 30, so 
n v/ n N ~ 10 11 . An order-of-magnitude estimate for the annihilation cross-section for the process 
Nv —> x a gives ( av } ~ y 2 /m 2 , where m ~ 10 — 100 GeV is a characteristic energy scale and 
y 2 < 10 -16 . Noting that the sneutrino self-annihilation cross-section at freeze-out is ~ 10~ 9 GeV~ 2 
(as required to obtain the right DM abundance), the ratio of the annihilation rate for N with v 
over that for N with N is ~ y 2 10 2 °/m 2 < 10 4 /m 2 . This shows that the process Nv —> x a nray be 
significant only close to the upper bound on y N determined above, but negligible for smaller values 
(i.e. y N < 10~ 9 — 10~ 10 ). However, if this reaction is below threshold, which is typically true given 
the metastability condition (5.15) and the comparable masses of the N and the x i n the region 
where the N freezes-out with the required abundance, the thermally-averaged cross-section will be 
much more suppressed. We do not seek to determine this bound on y N more precisely as it is clear 
that, at worst, this is only an order of magnitude stronger than that derived above and there is 
always available parameter space where y N can be made sufficiently small. 

6 Conclusions 

An excess in the 7 -ray spectrum originating from the galactic centre has been observed in the Ferrni- 
LAT data [5-10], and was recently confirmed by the Fermi-LAT collaboration [13]. Notwithstanding 
a possible astrophysical origin for this excess, it is well-explained by DM annihilating into 66 - or 
Higgs-pairs. Supersymmetric models provide a unifying framework to study the implications of 
DM annihilating via these channels since they contain a DM candidate and address many issues 
that arise beyond the SM. 

A suitable DM candidate in supersymmetric models is the neutralino. Since the annihilation of 
neutralinos into Higgs pairs is p -wave suppressed, we have considered the alternative channel of a 
Higgs and a Higgs-sector pseudoscalar. The latter predominantly decays into 66-pairs and thereby 
behaves similarly to the Higgs. We have performed a fit to the excess for this channel and found it 
to be just as good or even better than the 66- and Higgs-channel for a relatively light pseudoscalar in 
the final state, with mass below about 120 GeV. The quality of the fit prefers that the annihilation 
occurs as close to threshold as possible, but is typically acceptable to within 40 GeV. The best-fit 
cross-section lies within 0.4 — 30 x 10 -26 cm 3 /s, where the large range is mainly determined by 
uncertainties in the DM halo profile. 

In the MSSM the required light pseudoscalar is in strong tension with collider and flavour con¬ 
straints. We have therefore considered the NMSSM, which has an additional singlet superfield. 
Mixing of the singlet pseudoscalar with the MSSM pseudoscalar then eases the experimental con- 
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straints. In addition, in the regime of a large Higgs-singlet coupling A, the naturalness of the model is 
considerably improved. Direct-detection experiments limit the Higgsino fraction of the neutralino. 
Focusing instead on bino- and singlino-dominated neutralinos, we have performed a numerical 
scan over the NMSSM parameter space and found interesting regions that explain the excess via 
the Higgs-pseudoscalar channel. The points that we found do not rely on resonant enhancement 
of the annihilation cross-section, which was typically close to the canonical WIMP cross-section 
~ 2 x 10 _ 26 cm 3 /s. However, the masses m a and my of the pseudoscalar and the neutralino have to 
be chosen such that the annihilation is close to threshold (2my ~ m a + mh)- Furthermore, to avoid 
large deviations in the Higgs couplings from SM values, we found it necessary to have an accidental 
cancellation in the off-diagonal element of the Higgs mass matrix. Interestingly, the parameter 
space for the bino-like neutralino will soon be tested by direct-detection experiments, while the 
singlino-like neutralinos might escape direct detection even at the next-generation experiments. 

We have also performed a parameter scan for the 66 -channel, focusing on bino-dominated neu¬ 
tralinos. Due to direct-detection constraints on the coupling strength of the neutralino LSP to the 
Higgs sector, the annihilation cross section is typically not large enough to generate a sufficient 7 -ray 
flux for the excess. We have therefore relied on a light pseudoscalar in the s-channel with mass cho¬ 
sen such that the cross section is resonantly enhanced (m a — 2my; this mass relation may require a 
mild tuning). This results in a large increase of the cross section during freeze-out if 2my < m a ■ In 
this regime, all our points that could otherwise explain the excess therefore do not have a sufficient 
abundance to account for the observed amount of DM. In the opposite regime 2my > m a , we 
do find points which can explain both the excess and have a sufficient DM abundance. The vast 
majority of these points will be tested soon by the currently running direct-detection experiments 
and by improvements in constraining the invisible decay width of the Higgs. 

We have also studied an extension of the NMSSM with right-handed neutrinos that can explain 
the neutrino masses via the seesaw mechanism. The lightest right-handed sneutrino is then an 
interesting alternative DM candidate. It is straightforward to arrange for the sneutrino to annihilate 
dominantly into pseudoscalar-pairs. We have performed a fit for this channel and again found it to 
be just as good or even better than the 66 - and Higgs-channel for a relatively light pseudoscalar, 
with mass less than about 120 GeV. The best-fit cross-section is 0.3 — 30 x 10 ~ 26 cm 3 /s, where 
the large range is again dominated by uncertainties in the DM halo profile. Again it is preferred 
that the annihilation happens close to the threshold, typically within 60 GeV. We have performed 
a parameter scan of this extended NMSSM and found regions of parameter space where this new 
channel can be realised. It provides an interesting explanation of the excess as it can be easily 
obtained without relying on resonant enhancement of the annihilation cross-section. Finally, we 
have also considered a model where a right-handed sneutrino NLSP is metastable and decays after 
the freeze-out of the neutralino LSP. The relic abundance of the neutralino LSP can then be 
regenerated through these decays. The points for the 66 -channel with an underabundant neutralino 
LSP can thereby be resurrected. 

It is exciting that the 7 -ray excess from the galactic centre can be explained by DM annihilating 
into quark or scalar pairs and that this is straightforward to implement in the NMSSM with 
improved electroweak naturalness. Parts of the viable parameter space that we have found will soon 
be explored in direct-detection experiments. In addition, the 7 -ray line from scalar or pseudoscalar 
decays into photons may possibly be detected in the future. This would be an astonishing signal 
for supersymmetry from the sky. 


38 



Acknowledgements 

We thank Iason Baldes, Brian Batell, Nick Rodd, and Alfredo Urbano for useful discussions. This 
work was supported in part by the Australian Research Council. The work of TG was supported by 
the Department of Energy grant DE-SC0011842 at the University of Minnesota. We acknowledge 
the use of matplotlib [107], ipython [108] and Julia [109; 110]. 

A (S)Neutrino couplings 

The relevant couplings in the right-handed neutrino sector are given by 


C N R x k N ~ 

— V 2 \n A4 5 

(A.l) 

C NiX k N = 

—1V2 Atv A4.5 

(A.2) 

CaNN = 

— 1 \/ 2 \n Vis 

(A.3) 

ChNN = 

-V 2 Xn S13 

(A.4) 

C hN R N R = 

-V 2 X n (5i 3 (A AiV + 4 v s X n + 2 kv s ) - Xv (sin/35 n + cos/?5i 2 )) 

(A.5) 

C hN I N I = 

V 2 X N (Si 3 (A\ n - 4 v s X n + 2 nv s ) - Xv (sin/35n + cos/35i 2 )) 

(A.6) 

C hhN R N R = 

—2 A tv ((2A tv + k )Si 3 ~ A5i 2 5n) 

(A.7) 

C hhN I N I = 

—2 Xn ((2Aa t — K)Sf 3 + A5i 2 5n) 

(A.8) 

C aaN R N R 

—2 Xn ((2Aat — n)Vl 3 + A V12V11) 

(A.9) 

C aaN I N I = 

—2 A n ((2 Atv + k)'Pi3 ~ A V12V11) 

(A.10) 

C aN I N R = 

-V2X n ((A Ajv - 2 kv s )Vi 3 + A v (V12 cos /3 + Vu sin /3)). 

(A.ll) 


Our notation follows [73] and [35]. For NMSSM-specific couplings, see [73]. 
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